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This  report  describes  calculations  of  ionization  rales  and  electric  currents 
in  air  and  ground  resulting  from  neutron  sources.  In  particular,  this  report  des¬ 
cribes  the  numerical  methods,  nuclear  data,  and  computer  programs  employed 
and  presents  results  for  a  typical  thermonuclear  neutron  spectrum.  A  discussion 
of  fitting  of  earlier  electromagnetic  pulse  (EMP)  source  calculations  is  also 
presented. 

The  latest  ENDF/B  neutron  transport  and  gamma  ray  production  data 
for  air  and  ground  are  used  and  sensitivity  of  EMP  calculations  to  cross  section 
changes  have  been  estimated.  A  new  source  energy  band  version  of  the  Monte 
Carlo  code  05RNIES  has  been  developed  which  makes  it  possible  to  compute 
EMP  sources  by  folding  in  various  desired  neutron  spectra.  New  analytical 
approaches  have  also  been  introduced  to  accurately  treat  neutron  capture  events 
in  air  to  times  as  late  as  one  second.  Analytical  gamma-ray  transport  calcula¬ 
tions  away  from  the  air-ground  interface  are  also  used.  Various  neutron  and 
gamma-ray  contributions  to  ionization  rates  and  electric  currents  have  been 
tabulated  as  functions  of  fast  and  thermal  reactions  to  assist  in  fitting  of  these 
quantities. 

Calculations  are  presented  for  three  heights  of  burst,  0  meters,  2C0 
meters,  and  500  meters.  The  early  time  mesh,  below  0. 1  Msec,  is  more 
detailed  for  the  HoB  =  0  m  calculation  than  in  previous  calculations. 
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This  report  describes  calculations  of  ionization  rates  and  elec¬ 
tric  currents  in  air  and  ground  resulting  from  neutron  sources.  In 
particular,  this  report  describes  the  numerical  methods,  nuclear  data, 
and  computer  programs  employed  and  presents  results  for  a  typical  ther¬ 
monuclear  neutron  spectrum.  A  discussion  of  fitting  of  earlier  electromag¬ 
netic  pulse  (EMP)  source  calculations  is  also  presented. 

The  latest  ENDF/B  neutron  transport  and  gamma  ray  production 
data  for  air  and  ground  are  used  and  sensitivity  of  EMP  calculations  to 
cross  section  changes  have  been  estimated.  A  new  source  energy  band 
version  of  the  Monte  Carlo  code  05RNIES  has  been  developed  which  makes 
it  possible  to  compute  EMP  sources  by  folding  in  various  desired  neutron 
spectra.  New  analytical  approaches  have  also  been  introduced  to  accu¬ 
rately  treat  neutron  capture  events  in  air  to  times  as  late  as  one  second. 
Analytical  gamma- ray  transport  calculations  away  from  the  air-ground 
interface  are  also  used.  Various  neutron  and  gamma-ray  contributions 
to  ionization  rates  and  electric  currents  have  been  tabulated  as  functions 
of  fast  and  thermal  reactions  to  assist  in  fitting  of  these  quantities. 

Calculations  are  presented  for  three  heights  of  burst,  0  meters, 
200  meters,  and  500  meters.  The  early  time  mesh,  below  0. 1  ^sec,  is 
more  detailed  for  the  HoB  =  0  rn  calculation  than  in  previous  calculations. 
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1.  INTRODUCTION 


Calculations  are  presented  of  electromagnetic  pulse  (EMP) 
sources  for  a  typical  thermonuclear  neutron  source  spectrum.  The 
ionization  rates  and  Compton  currents  resulting  from  the  neutron 
source  are  required  for  the  prediction  of  free  field  EMP  environments. 
The  calculations  are  performed  for  a  point  source  of  neutrons  in  an  air- 
over-ground  configuration.  Three  low  level  heights  of  burst  are  con¬ 
sidered;  0,  200,  and  500  meters.  A  time- dependent  tabulation  of  ioni¬ 
zation  rates  and  Compton  currents  is  made  in  two  dimensions  in  both 
air  and  ground. 

New  analytical  methods  have  been  developed  to  take  into  con¬ 
sideration  the  very  late-time  neutron  capture  reactions  in  air.  These 
events,  which  can  introduce  sources  at  times  as  late  as  one  second, 
have  been  included  in  our  calculations  for  the  first  time.  A  description 
of  the  new  version  of  the  Monte  Carlo  code  05RNIES  ,  which  incorpor¬ 
ates  these  air  capture  reactions,  is  included  in  a  later  section.  A  new 
energy  band  treatment  enables  the  user  to  efficiently  fold  in  an  arbitrary 
neutron  spectrum  once  the  initial  production  runs  are  completed. 


(1_3) 

The  calculations  presented  here  are  similar  to  previous  results.'  ’ 
A  complete  treatment  of  neutron  transport  in  both  air  and  ground  is  in¬ 
cluded.  The  inclusion  of  late-time  air  capture  reactions  represents  a 
significant  improvement  over  earlier  treatments  which  neglected  these 
events.  Thermal  neutrons  are  also  followed  back  and  forth  from  ground 
to  air  in  the  new  calculations,  as  opposed  to  previous  results  which  ap¬ 
proximated  the  air  by  vacuum.  Comparisons  made  between  the  new  late- 
time  program  and  the  previous  version  of  05RNIES  indicate  that  the  new 
version  yields  substantially  higher  ionization  rates  than  the  conventional 
program  in  the  time  interval  between  1  millisecond  and  1  second.  In  this 

1 


time  span  air  capture  events,  which  are  neglected  in  the  earlier  version 
of  05RNIES,  become  important  and  enhance  the  late-time  ionization 
rates  as  computed  in  the  new  code.  Furthermore,  ionization  rates  and 
currents  in  the  ground,  but  very  close  to  the  air -ground  interface,  also 
appear  to  be  enhanced  relative  to  the  conventional  treatment.  This  is 
due  in  part  to  the  fact  that  the  conventional  version  of  the  code  does  not 
allow  thermal  neutrons  escaping  from  the  ground  into  the  air  to  return 
to  the  ground,  whereas  the  late -time  models  follow  the  neutrons  from 
one  medium  into  another. 

Another  important  change  has  been  made  in  the  05RNIES  code. 
Previously,  only  the  total  ionization  rate  and  Compton  current  were  edited 
as  output.  In  the  new  version,  contributions  are  separated  according  to 
their  origin.  Thus,  for  current  rate  we  tabulate,  in  addition  to  the  total, 
the  separate  contributions  from  thermal  neutron  capture  gamma  rays  and 
from  fast  neutron  capture  and  inelastic  gamma  rays.  For  ionization  rate, 
the  contribution  from  neutron  collisions  is  also  tabulated  separately  as 
fast  and  thermal  neutron  induced  ionization.  These  separate  contributions 
will  give  greater  understanding  of  which  reaction  types  are  most  important 
at  different  times.  This  additional  knowledge  is  most  useful  when  selecting 
biasing  techniques  and  will  make  it  much  easier  to  fit  the  results  because 
the  behavior  of  the  separate  components  is  easier  to  describe  with  simple 
functions  than  is  the  sum  of  the  components. 

The  results  presented  here  incorporate  the  latest  ENDF/B  neutron 
transport  and  gamma  ray  production  cross  sections  for  air  and  ground. 
Sensitivity  of  EMP  calculations  to  changes  in  cross  sections  has  also  been 
estimated. 

The  calculations  are  carried  out  to  local  times  of  one  second  and 
to  ranges  of  five  kilometers.  However,  the  time -dependent  data  at  ranges 
beyond  about  three  kilometers  are  somewhat  erratic  with  large  statistical 
uncertainties. 
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2.  DISCUSSION  AND  RESULTS  FOR  A  TYPICAL 
THERMONUCLEAR  NEUTRON  SPECTRUM 

EMP  source  calculations  are  presented  in  this  section  for 
three  low  altitude  nuclear  explosions.  The  three  heights  of  burst  (HoB's) 
selected  were  at  0,  200,  and  500  meter  altitudes.  The  neutron  source 
spectrum  chosen  for  these  calculations  was  a  typical  thermonuclear 
spectrum  which  has  been  used  in  other  transport  calculations. '  The 
typical  thermonuclear  spectrum  is  shown  in  Table  1  in  a  seventeen  group 
structure.  The  actual  calculations  were  performed  using  thirteen  energy 
bands  as  will  be  discussed  in  Section  4.2. 

We  have  produced  a  total  of  more  than  5,000  ionization  rate 
and  current  plots  on  microfilm  for  the  thermonuclear  source.  It  is 
apparent  that  we  can  only  present  a  very  small  fraction  of  these  results 
in  this  report.  However,  we  do  present  a  representative  sampling  of 
these  data  and  have  stored  the  complete  plotted  and  tabulated  results  on 
microfilm  and  tape  where  they  may  be  processed  for  use  in  EMP  field 
codes.  Furthermore,  the  basic  production  runs  are  available  on  tape  for 
each  height  of  burst  so  that  other  neutron  spectra  may  easily  be  folded  in 
to  obtain  various  desired  EMP  drivers. 

Ionization  rates  and  Compton  currents  were  obtained  in  air  and 
ground  and  extend  to  ranges  of  5  kilometers  in  air.  Biasing  in  the  form 
of  path  length  stretching  of  neutrons  and  gammas  is  mandatory  to  obtain 
results  at  these  distances  with  reasonable  amounts  of  computer  time. 
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Table  1 


ENERGY  DISTRIBUTION  FOR 
THERMONUCLEAR  SOURCE 


Energy 

(MeV) 

Number 
Fraction 
in  Grot?) 

12.2  - 

15.0 

7.  06(-2)J 

10.0  - 

12.2 

2.  56(-2) 

8.19  - 

10.0 

1.41  (-2) 

6.36  - 

8.19 

1.47(-j) 

4.97  - 

6.36 

1 . 80(-2) 

4.07  - 

4.97 

1 . 70(-2) 

3.01  - 

4.07 

2.60(-2) 

2.46  - 

3.01 

1 . 90(-2) 

2.35  - 

2.46 

5.  00(-3) 

1.83  - 

2.35 

2.  80(-2) 

1.11  - 

1.83 

6. 20(-2) 

0.55  - 

1.11 

8.  50(-2) 

0.11  - 

0.55 

1 . 02  (-  1) 

3. 35 (-3)  - 

0.110 

3.  65(- 1) 

5. 83 (-4)  - 

3.  35(-3) 

1 . 22(-  1) 

1 . 01  (-4)  - 

5.  83(-4) 

2.  40(-2) 

2.  90(-5)  - 

1 . 01 (-4) 

2.  00(-3) 

<2.  90(-5) 

0 

^ead  as  7. 06  x  10~2 


u 

Even  with  biasing,  the  spatial  results  beyond  3  kilometers  in  air,  and 
also  at  great  depths  in  the  ground,  become  sporadic.  However,  it  is 
felt  that  most  of  the  spatial  positions  of  interest  in  EMP  studies  have 
been  obtained  with  reasonable  statistical  accuracy. 

The  calculations  described  here  were  acquired  with  a  new 
late-time  version  of  the  Monte  Carlo  code  05RNIES,  which  is  discussed 
in  detail  in  Section  4.  The  latest  ENDF/B  neutron  and  gamma  ray  pro¬ 
duction  data  were  used  and  the  sensitivity  of  calculations  to  cross  section 
changes  is  considered  in  Section  4.  5.  The  accurate  treatment  of  very 
late-time  air  capture  reactions  represents  a  significant  improvement  in 
the  calculational  models  used  in  generating  EMP  sources.  The  new  late¬ 
time  results  are  generally  higher  than  previous  results  in  air  in  the  time 
range  from  1  millisecond  to  1  second,  where  air  capture  events  become 
important.  The  accurate  tracking  of  thermal  neutrons  back  and  forth 
from  ground  to  air  in  the  new  code  also  generally  results  in  larger  EMP 
sources  at  positions  in  the  ground  but  near  the  air -ground  interface. 

The  05RNIES  code  tabulates  results  in  spherical  coordinates 
above  the  source,  and  in  cylindrical  coordinates  below  it.  This  coordi¬ 
nate  system  is  selected  for  consistency  with  the  EMP  field  codes  which 
use  the  05RNIES  results  as  input. 

In  each  of  the  three  HoB  calculations  discussed  here,  a  total 
of  195,000  neutron  histories  were  followed.  These  neutrons,  through 
various  inelastic  scattering,  capture,  charged  particle,  and  other  reac¬ 
tions,  produce  gamma  rays  which  are  also  followed.  For  each  HoB, 
production  runs  were  performed  with  unit  neutron  sources  in  each  of 
thirteen  energy  bands.  The  desired  neutron  source  spectra  were  then 
folded  in  using  the  auxiliary  folding  and  plotting  code  05NPT.  The  pro¬ 
duction  rims  required  approximately  22  hours,  17  hours,  and  15  hours 
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of  computer  time  on  the  CDC  6600  for  the  HoB  =  0  m,  HoB  -  200  m,  and 
HoB  =  500  m  runs,  respectively.  The  higher  altitude  bursts  are  per¬ 
formed  more  rapidly  since  the  new  faster  analytic  routines  in  05RNIES 
can  be  employed  more  frequently  at  the  higher  elevations.  The  auxiliary 
folding  code  05NPT  required  approximately  30  minutes  of  CDC  6600  time 
to  fold  in  an  arbitrary  neutron  spectrum  for  each  HoB. 

The  bin  structures  used  for  scoring  for  the  three  heights  of 
burst  are  indicated  in  Tables  2-4.  The  results  are  obtained  to  dis¬ 
tances  of  5  kilometers  and  to  times  of  1  second.  For  the  HoB  =  0  m 
calculations,  there  are  39  x  19  x  14  =  10,374  individual  scoring  bins. 

Similarly,  for  the  HoB  =  200  m  and  HoB  =  500  m  problems  we  score 
time  dependent  quantities  in  9120  and  9690  bins,  respectively.  It  is 
apparent  that  with  the  large  number  of  bins  in  each  problem,  good  sta¬ 
tistical  results  could  not  be  obtained  at  each  position.  However,  it  is 
felt  that  all  of  the  important  spatial  positions  are  adequately  represented 
in  the  results  described  here. 

Time  integrated  ionization  data  are  considered  in  the  first 
twelve  figures.  Following  these,  time  dependent  ionization  rate  and 
Compton  current  calculations  are  presented.  All  of  the  figures  are 
normalized  to  one  source  neutron.  In  some  of  the  figures,  contributions 
due  to  thermal  and  fast  neutron  reactions  are  separately  plotted.  Be¬ 
cause  of  space  limitations,  the  statistical  error  bars  which  depict  stan¬ 
dard  deviations  from  the  mean  values  are  only  shown  for  the  total 
quantities.  The  separation  into  fast  and  thermal  components  was  made 
to  assist  in  fitting  of  these  data  and  to  further  understanding  of  the  im¬ 
portance  of  various  reaction  mechanisms  in  different  time  regimes. 

In  Fig.  2. 1  the  ionization  vs. radius  is  shown  for  an  angular  bin 
having  polar  angle  limits  of  cos  8  =  0. 0  to  cos  8  =  0. 05  for  a  source  on 
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the  ground.  The  ionization  is  due  to  both  neutrons  and  gamma  rays  and 
the  solid  histogram  depicts  tte  total  dose.  The  dashed  histogram  repre¬ 
sents  contributions  to  the  ionization  from  fast  neutron  ionization  processes 
and  from  gamma  rays  which  are  produced  from  fast  neutron  reactions. 

The  dashed-dotted  histogram  similarly  indicates  ionization  produced  by 
thermal  neutron  energy  deposition  and  by  photons  which  are  produced 
by  thermal  neutron  reactions. 

Figure  2. 2  shows  the  ionization  vs.  radius  for  the  200  m  HoB 
for  a  height  bin  extending  from  the  air-ground  interface  to  an  elevation 
of  50  meters.  These  bins  represent  approximately  the  same  spatial 
positions  as  the  angular  bins  of  Fig.  2. 1,  namely  spatial  regions  close 
to  the  air -ground  interface.  Comparing  the  0  and  200  meter  HoB 
results  along  the  air-ground  interface  indicates  that  there  are  differences 
in  the  time  integrated  ionization  rate  results.  At  close-in  ranges  out 
to  a  few  hundred  meters, the  0  meter  results  are  substantially  higher 
while  at  greater  distances  the  200  meter  HoB  results  appear  to  be 
generally  higher. 

In  Fig.  2.3  the  time  integrated  ionization  rate  vs.  radial  position 
is  shown  for  the  500  meter  HoB  for  the  elevation  bin  extending  from  the 
ground  to  an  altitude  of  50  meters.  The  500  meter  HoB  results  are 
much  lower  than  the  200  meter  HoB  calculations  for  the  first  few  radial 
bins  extending  to  several  hundred  meters.  The  500  and  200  meter  HoB 
results  are  then  fairly  close  to  one  another,  with  the  200  meter  results 
generally  higher,  out  to  about  2  kilometers,  beyond  which  the  500  meter 
calculations  surpass  the  lower  HoB  data. 

A  comparison  of  the  0  and  200  meter  HoB  time  integrated 
ionization  rates  vs.  radius  for  the  angular  bin  cos  6=  0.5  -1.0  is 
shown  in  Fig.  2.  4.  The  values  for  the  HoB  =  0  source  are  generally 
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Figure  2.2 


Radial  dependence  of  ionization  due  to  thermo¬ 
nuclear  source  for  HoB  =  200  m,  altitude  = 
0-50  meters. 


higher  than  the  HoB  =  200  m  results  for  all  ranges.  A  similar  comparison 
was  made  between  the  200  meter  HoB  and  the  500  meter  HoB  for  the 

cos  6  =  0.  5  -  1. 0  bins.  For  these  two  cases  virtually  identical  results 

* 

were  obtained  for  all  ranges  which  implies  that  for  HoB's  above  200 
meters,  the  ground  plays  little  role  in  determining  the  dose  deposited  in 
air  at  small  polar  angles. 

The  reason  for  the  differences  between  the  0  and  200  meter  HoB 
results  shown  in  Fig.  2. 4  can  be  explained  with  the  aid  of  Fig.  2. 5-2. 11. 
From  Fig.  2.5  it  is  clear  that  secondary  photons  which  are  produced 
in  the  ground  by  a  HoB  =  0  meter  source  must  traverse  more  ground 
to  reach  an  observer  at  b  than  an  observer  at  a.  Therefore,  the  energy 
deposited  by  gamma  rays  at  small  polar  angles  in  air  should  be  larger 
than  the  energy  deposited  at  larger  angles.  A  similar  argument  can  be 
presented  for  the  source  neutrons  from  a  ground  burst.  A  large  fraction 
of  the  source  neutrons  produce  secondary  neutrons  in  the  ground  through 
inelastic  scattering  reactions.  The  energy  deposited  by  these  neutrons 
is  again  larger  for  small  polar  angles  because  there  is  less  ground  to 
traverse  than  for  the  large  polar  angles. 

The  polar  angle  dependence  of  the  ionization  for  the  ground 
burst  is  shown  in  Figs.  2.  6-2.  8.  In  Fig.  2.  6  the  total  neutron  plus 
photon  ionization  as  a  function  of  the  polar  angle  is  presented  for  a  radial 
bin  extending  from  450  to  550  meters.  The  ionization  due  to  both  fast 
and  thermal  reactions  is  seen  to  increase  with  decreasing  polar  angle 
as  would  be  expected  from  the  previous  geometrical  arguments  of 
Fig.  2.5.  Figures  2.  7  and  2.  8  illustrate  the  ionization  versus  polar 
angle  due  to  neutrons  and  gamma  rays,  respectively.  For  both  types  of 
particles,  it  is  clear  that  the  ionization  increases  with  decreasing  polar 
angle.  Neutron  and  photon  ionization  produced  by  both  fast  and  thermal 
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Figure  2. 5. 


Geometrical  explanation  for  polar  angle 
dependence. 
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Figure  2  7.  Polar  angle  variation  of  neutron  ionization  at 
radius  450  -  550  meters,  due  to  thermonuclear 
source  on  ground. 
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Figure  2. 8.  Polar  angle  variation  of  gamma  ray  ionization  at 
radius  450  -  550  meters,  due  to  thermonuclear 
source  on  ground. 


neutron  reactions  demonstrate  this  angular  dependence.  Although  these 
angular  results  are  shown  here  only  for  one  radial  bin,  the  effect  can  be 
shown  to  be  a  general  phenomenon  at  other  ranges  for  a  ground  burst.  For 
higher  HoB's,  the  polar  angle  dependence  is  much  flatter.  In  Figs.  2.9- 
2. 11,  the  polar  angle  dependence  of  the  ionization  at  the  450-550  meter 
radial  bin  is  shown  for  the  200  meter  HoB.  The  total  ionization,  neutron 
ionization,  and  photon  ionization  results  are  all  flatter  for  the  200  meter 
HoB  as  a  function  of  polar  angle  than  the  corresponding  ground  burst  calcu¬ 
lations.  The  polar  angle  dependence  is  much  less  pronounced  for  both  the 
200  m  and  500  m  HoB  sources  than  for  the  ground  burst  because  for  the 
higher  altitude  calculations  the  secondary  photons  and  neutrons  produced 
in  the  ground  are  distributed  over  a  wider  area  and  furthermore,  the 
ground  plays  a  much  smaller  role  in  altering  directions  for  source 
particles  because  of  the  greater  distances  involved.  Therefore,  for 
small  polar  angles,  larger  ionization  values  due  to  the  ground  burst, 
relative  to  the  200  m  HoB,  can  be  anticipated  as  illustrated  in  Fig.  2.4. 

Figure  2. 12  illustrates  the  effect  of  the  height  of  the  source  upon 
ionization  in  the  ground  for  the  radial  interval  from  450-550  meters. 

The  HoB  =  200  m  results  are  consistently  higher  than  the  ground  burst 
calculations.  In  turn,  the  HoB  =  0  m  results  are  slightly  higher  than 
the  higher  altitude  HoB  =  500  m  data.  The  HoB  =  200  m  calculations  exceed 
the  ground  burst  mainly  due  to  the  fact  that  uncollided  neutrons  can  reach 
the  ground,  at  this  radial  distance,  from  the  200  m  burst  but  not  from  the 
ground  burst.  Therefore,  the  HoB  =  0  m  source  neutrons  are  degraded 
in  energy  to  a  greater  extent  than  the  200  m  particles,  and  the  ionization 
in  the  ground  is  therefore  substantially  reduced  relative  to  the  HoB  =  200  m 
calculations.  On  the  other  hand,  as  even  higher  HoB  sources  are  con¬ 
sidered,  the  probability  of  energetic  neutrons  and  gamma  rays  reaching 
the  ground  at  the  450-550  m  radius  is  reduced.  This  is  demonstrated  by 
the  fact  that  the  HoB  =  500  m  results  are  slightly  lower  at  nearly  all  depths 
than  the  ground  burst  results,  and  much  lower  than  the  HoB  =  200  m 
calculations. 
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Figure  2. 10.  Polar  angle  variation  of  neutron  ionization  at 
radius  450-550  meters,  due  to  thermonuclear 
source  for  200  m  HoB. 
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Figure  2. 12.  Ionization  vs.  depth  at  radius  450  -  550  meters, 
due  to  thermonuclear  source  for  HoB  =  0,  HoB  = 
200,  and  HoB  =  500  meters. 


The  next  series  of  figures  is  concerned  with  various  time 
dependent  ionization  rate  and  Compton  current  results  at  radial  positions 
of  450  -  550  meters  for  the  three  heights  of  burst.  In  Fig.  2.13  the 
ionization  rates  in  air  just  above  the  air -ground  interface  are  displayed 
for  the  HoB  =  0  m  and  HoB  =  200  m  calculations.  The  HoB  =  200  m 
results  are  for  a  cylindrical  bin  extending  from  the  ground  to  an 
elevation  of  50  meters.  The  HoB  =  0  m  results  are  for  a  spherical 
bin  extending  from  cos  0  =  0. 0  to  cos  0  =  0. 05.  Although  the  two 
bins  are  of  different  shapes,  they  both  refer  to  approximately  the  same 
spatial  positions.  The  results  presented  in  Fig.  2. 13  indicate  that  there 
is  very  little  difference  between  the  two  HoB  calculations  for  the  spatial 
positions  considered  here.  The  HoB  *  0  m  results  are  higher  in  the 
very  early  time  intervals  since  neutrons  are  immediately  converted  to 
inelastic  gammas  for  the  ground  source.  Therefore,  the  HoB  =  0  m 
calculations  display  an  initial  prompt  spike  which  is  absent  in  the  higher 
altitude  calculations.  Both  sources  show  another  spike  at  approximately 
8  fisec  local  time  which  signals  the  arrival  of  14  MeV  neutrons.  After 
the  arrival  of  the  neutron  wave  front,  ionization  is  caused  by  both 
photons  and  neutrons,  whereas  at  earlier  times  only  secondary  gamma 
rays  were  contributing  to  ionization. 

In  Fig.  2. 14  a  comparison  between  HoB  =  200  m  and  HoB  =  500  m 
ionization  rates  is  made  for  the  elevation  bin  0-50  meters.  Here  the 
500  meter  results  are  lower  than  the  HoB  =  200  m  calculations  by  about 
a  factor  of  2  for  all  the  time  intervals.  The  arrival  of  the  14  MeV 
neutron  spike  is  also  delayed  for  the  HoB  =  500  m  calculations,  falling 
into  the  10  -  20  ftsec  local  time  bin,  since  these  neutrons  must  travel 
farther  before  reaching  the  observation  bin  than  the  lower  altitude  source 
neutrons. 


25 


Ionization  Rate  (MeV/m  -sec) 


10  “  10 


103  104 


Local  Time  (microseconds) 


Figure  2.13.  Ionization  rate  for  thermonuclear  source  at 
radius  =  450-550  meters,  for  HoB  =  0  m, 
cos  0  =  0-0. 05,  and  at  HoB  =  200  m, 
elevation  =  0  -  50  m. 
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Figure  2. 14.  Ionization  rate  at  radius  450-550  meters,  elevation 
=  0-50  m,  for  thermonuclear  source  at  HoB  =  200 
and  HoB  =  500  m. 


Figure  2.15  shows  the  ionization  rate  in  the  air  for  a  cos  6  = 

0. 5  -  1. 0  bin  due  to  a  ground  burst.  The  ionization  rate  is  due  to  combined 
neutron  and  gamma  ray  energy  deposition.  Figures  2.16  and  2.17  display 
the  ionization  rates  for  this  same  angular  bin  due  to  separate  neutron  and 
gamma  ray  deposition  processes.  On  all  three  plots,  the  ionization  rate 
is  further  separated  into  components  due  to  fast  and  thermal  neutron 
reactions.  By  comparing  the  three  figures,  we  can  determine  something 
about  the  effectiveness  of  the  various  deposition  mechanisms  in  certain 
time  intervals.  It  is  apparent  from  Figs.  2.16  and  2. 17  that  the  ionization 
is  solely  due  to  gamma  rays  initially.  A  prompt  initial  photon  spike  due  to 
inelastic  neutron  scattering  in  the  ground  is  displayed  in  the  first  time 
interval.  Photon  ionization  due  to  thermal  neutron  captures  does  not 
appear  until  approximately  4  nsec  local  time.  The  initial  fast  neutron 
energy  deposition  appears  at  about  5  nsec,  with  significant  deposition 
taking  place  in  the  next  time  bin,  from  7.  5  -  10  ^sec.  From  this  time  bin 
out  to  approximately  100  nsec,  the  neutron  ionization  is  larger  than  the 
photon  contribution.  The  gamma  ray  ionization  tnen  becomes  larger  than 

4 

the  neutron  component  in  the  time  range  from  100  to  10  nsec.  At 

4 

10  nsec  thermal  neutron  capture  in  aii  becomes  a  factor,  as  can  be  seen 

from  Fig.  2.16,  and  the  neutron  and  photon  components  of  the  ionization 

4  6 

rate  become  comparable  from  10  to  10  nsec. 

Figures  2.18-2.20  show  the  total  ionization  rate,  and  the  ioniz¬ 
ation  rates  due  to  the  separate  neutron  and  gamma  ray  processes  for 
a  source  at  an  altitude  of  500  m.  The  angular  bin  extending  from  cos  6  = 

0.  5  -  1.  0  with  radius  450  -  550  meters  is  again  used.  In  comparing  the 
HoB  =  500  m  results  with  the  ground  burst  calculations  of  Figs.  2.15- 
2. 17  we  can  draw  a  number  of  conclusions.  Neutron  ionization  rates  for 
the  two  sources  are  very  similar  for  all  time  bins  as  might  be  anticipated 
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Figure  2.16.  Ionization  rate  due  to  neutrons  at  radius  450  -  550 
meters,  cos  6  =  0. 5  -  1. 0,  due  to  thermonuclear 
source  at  HoB  =  0  meters. 
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Figure  2. 17.  Ionization  rate  due  to  gamma  rays  at  radius  450- 
550  meters,  cos  d  =  0.  5  -  1. 0,  due  to  thermo¬ 
nuclear  source  at  HoB  =  0  meters. 
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Figure  2.18.  Total  ionization  rate  at  radius  450  -  550  meters, 
cos  6  =  0. 5  - 1. 0,  due  to  thermonuclear  source  at 
HoB  =  500  meters. 
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Figure  2.19.  Ionization  rate  due  to  neutrons  at  radius  450 -  550 
meters,  cos  6  =  0. 5  -  1. 0,  due  to  thermonuclear 
source  at  HoB  =  500  meters. 
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Figure  2. 20.  Ionization  rate  due  to  gamma  rays  at  radius 
450  -  550  meters,  cos  8  =  0. 5  -  1. 0,  due  to 
thermonuclear  source  at  HoB  =  500  meters. 
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since  the  spatial  bin  considered  is  in  the  air  and  distant  from  the  air- 
ground  interface.  The  higher  altitude  calculations  again  do  not  display 
the  large  prompt  spike  shown  by  the  HoB  =  0  m  calculations.  Further¬ 
more,  the  ground  burst  shows  substantially  higher  ionization  out  to  about 
0.4  (i  sec  where  the  two  calculations  are  approximately  equal.  This  is 
expected  since  the  gamma  rays  from  inelastic  scattering  in  the  ground 
are  only  possible  from  the  HoB  =  0  m  source  in  this  time  region.  The 
HoB  =  500  m  fast  gamma  rays  are  due  solely  to  inelastic  scattering  in 
the  air  out  to  about  10  m  sec  where  ground  inelastic  gammas  become  a 
possibility.  The  ionization  due  to  photons  is  approximately  the  same 

4 

for  the  two  HoB's  from  .0. 4  m sec  to  3  nsec.  From  3  psec  to  about  10 

Msec  the  HoB  =  0  m  results  are  about  an  order  of  magnitude  larger 

4 

than  the  500  m  calculations.  Beyond  10  Msec  the  two  results  are 
comparable.  It  is  interesting  to  note  that  a  very  large  difference  in 
photon  ionization  exists  between  the  two  calculations  due  to  thermal 
photon  contributions.  The  ground  burst  thermal  photon  ionization  is  about 
three  orders  of  magnitude  larger  than  the  HoB  =  500  m  calculation.  The 
thermal  gamma  contributions  also  appear  at  earlier  times,  4  Msec  for 
the  ground  burst  versus  30  psec  for  the  higher  source.  Differences  in 
thermal  gamma  ionization  rates  between  the  two  sources  appear  to  be 
primarily  responsible  for  the  large  deviations  displayed  in  the  total 
gamma  ionization  rates. 

In  Fig.  2. 21  the  radial  Compton  current  due  to  the  ground  burst  is 
shown  for  the  angular  bin  extending  from  cos  6  =  0. 0  to  cos  6  =  0. 05 
and  the  range  from  450  -  550  m.  Figure  2. 22  displays  the  theta  current 
for  the  same  angular  bin.  In  both  figures  negative  currents  are  indicated 
by  the  letter  "N”.  The  sign  convention  used  on  all  figures  in  this  report 
is  that  the  electric  current  flows  in  the  direction  opposite  to  the  electron 
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Figure  2. 22.  Theta  current  at  radius  450 -  550  meters, 

cos  6  =  0-0. 05, due  to  thermonuclear  source 
on  ground. 
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flow.  By  including  both  air  and  ground  capture  gamma  production  in  our 
calculations,  electric  currents  are  produced  to  much  longer  times  than 
in  previous  studies  v  .  The  radial  currents  in  Fig.  2.21  are  essentially 
all  negative  as  would  be  anticipated  indicating  that  the  photon  current 
is  in  an  outward  direction  at  all  times.  On  the  other  hand,  the  theta 
current  is  seen  to  undergo  two  distinct  sign  reversals.  At  early  times 
to  approximately  100  nsec,  the  photon  current  is  predominantly  into  the 
ground  resulting  in  a  negative  theta  current.  For  times  extending  to 

4 

about  3  x  10  nsec,  gamma  rays  are  produced  primarily  in  the  ground 
and  the  transport  of  these  photons  into  the  air  produces  positive  theta 
currents.  At  very  late  times,  air  capture  gamma  rays  dominate,  and 
photon  currents  directed  toward  the  ground  again  result  in  negative  theta 
currents.  This  very  late  time  sign  reversal  of  the  theta  current  is  a 
new  phenomenon  which  was  not  observed  in  last  year’s  calculations 
since  air  capture  events  were  not  considered  at  that  time. 

Figures  2.23  and  2. 24  show  the  radial  currents  for  the  angular 
bins  cos  6  =  0. 0  -  0. 1  and  radial  distances  450  -  550  m  due  to  the 
HoB  =  200  m  and  HoB  =  500  m  sources,  respectively.  The  currents  for 
the  two  sources  are  in  reasonable  agreement  for  nearly  all  time  bins.  Further¬ 
more  comparing  the  200  m  HoB  with  the  ground  burst  of  Fig.  2. 21  again 
reveals  good  agreement  for  most  time  intervals.  The  ground  burst  of 
Fig.  2.21  again  displays  the  initial  prompt  gamma  spike  due  to  inelastic 
collisions  in  the  ground,  whereas  the  two  higher  altitude  bursts  are 
smoother  and  of  lower  magnitude  at  very  early  times.  Comparison  of  the 
theta  currents  for  the  two  higher  HoB’s  with  the  ground  burst  shows 
much  larger  discrepancies  which  may  be  partly  due  to  larger  statistical 
uncertainties.  However,  the  theta  currents  for  the  higher  bursts  also 
display  sign  changes  similar  to  those  of  Fig.  2.22,  although  the  changes 
take  place  in  different  time  bins  for  each  source. 
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Figure  2.23.  Radial  current  at  radius  450 -  550  meters, 

cos  6  =  0  -  0. 1,  due  to  thermonuclear  source  for 
HoB  =  200  meters. 


Local  Time  (microseconds) 


Radial  current  at  radius  450  -  550  meters, 
cos  0  =  0  -  0. 1,  due  to  thermonuclear  source  for 
HoB  =  500  meters . 


In  Fig.  2. 25  a  comparison  is  made  of  the  ionization  rate  for  the 

cos  8  =  0. 5  - 1. 0  bin  between  the  HoB  =  0  m  and  HoB  =  200  m  sources.  ' 

A  radial  range  of  450  -  550  m  is  again  used.  The  ionization  rate  at  early 

times  below  1  n sec  is  again  much  higher  for  the  ground  burst.  From 

1  //sec  to  100  //sec  the  two  sources  give  very  comparable  results,  while 

4 

between  100  and  10  /.sec  the  HoB  =  0  m  source  is  again  substantially 

4 

higher.  The  results  to  10  / isec  are  basically  consistent  with  conclusions 

(3  \  1 

drawn  in  comparing  the  two  sources  in  last  year's  effort.  v  '  However, 
the  deviations  between  the  two  sources  in  the  time  interval  100  - 10^  //sec 

appear  to  be  somewhat  larger  in  the  previous  work.  Furthermore,  j 

4 

beyond  10  //sec  the  two  sources  yield  approximately  equal  results.  It 

appears  that  the  inclusion  of  air  capture  events  in  the  present  calculations 

enhances  the  HoB  =  200  m  results  relative  to  the  ground  burst  calculations. 

We  have  also  compared  the  HoB  =  200  m  source  results  with  HoB  = 

500  m  calculations  for  this  angular  bin  and  found  close  agreement  between 

4  1 

the  two  sources  out  to  100  //sec.  From  100  // sec  to  3  x  10  //sec,  factor 
of  2  variations  are  observed  with  the  200  m  HoB  results  higher,  and 
good  agreement  is  again  seen  at  very  late  times.  In  general  the  agree¬ 
ment  between  the  HoB  =  200  m  and  HoB  =  500  m  sources  is  very  good 
for  this  spatial  bin. 

The  effects  of  the  neutron  source  altitude  on  ionization  rates  in 
the  ground  is  illustrated  in  the  next  four  figures.  In  Fig.  2. 26  a  com¬ 
parison  is  made  between  ionization  rates  due  to  a  ground  burst  and  due  to 
the  HoB  =  200  m  source  for  the  450  -  550  m  radial  bin  and  for  a  depth  of 
0-2  cm.  Beyond  the  first  few  fractions  of  a  microsecond,  the  HoB  = 

200  m  rates  appear  to  be  generally  higher  than  the  ground  source  results. 

A  comparison  between  HoB  =  200  m  and  HoB  =  500  m  ionization  rates 
for  the  same  ground  bin  is  presented  in  Fig.  2. 27.  Here  the  500  m  HoB 
is  consistently  lower  than  the  200  m  source  at  all  times.  For  this  bin, 
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Figure  2. 25.  Ionization  rate  at  radius  450  -  550  meters, 

cos  6  =  0. 5  - 1. 0,  due  to  thermonuclear  source 
at  HoB  =  0  and  HoB  =  200  meters. 


Ionization  Rate  (MeV/m  -sec) 


the  ionization  rate  due  to  the  HoB  =  200  m  source  exceeds  the  ground 
burst  results  since  neutrons  from  the  higher  altitude  source  can  reach 
the  observation  point  without  prior  attenuating  collisions  in  the  ground. 

The  ground  burst  source  neutrons  on  the  other  hand  are  moderated  to  a 
greater  extent  by  the  ground  at  positions  which  are  distant  from  the 
observation  bin.  Therefore,  there  is  less  neutron  energy  available  for 
deposition  in  this  bin.  The  probability  of  photon  energy  deposition  is 
also  reduced  because  the  photons  are  born  at  positions  which  have  smaller 
solid  angle  probabilities  of  eventually  reaching  the  observation  point. 

The  HoB  =  500  m  ionization  rates  are  generally  lower  than  the  200  m 
source  results  because  the  higher  altitude  source  neutrons  are  attenuated 
in  the  air  to  a  greater  extent  than  the  HoB  =  200  m  neutrons.  Therefore, 
the  500  m  HoB  neutron  energy  available  for  depos  '  ion  in  the  0  -  2  cm  bin 
is  less  than  the  neutron  energy  available  for  the  lower  altitude  source. 
Furthermore,  the  gamma  ray  energy  deposited  by  the  500  m  HoB  source 
is  lower  than  the  HoB  =  200  m  results  due  to  greater  attenuation  of 
gamma  rays  which  are  produced  at  the  higher  altitudes. 

In  Figs.  2. 28  and  2. 29  a  comparison  of  the  three  source  altitudes 
is  made  for  the  same  range  of  450  -  550  m,  but  for  a  deeper  ground  depth 
of  90  -  120  cm.  Figure  2.  28  compares  the  ground  burst  and  the 
HoB  =  200  m  source  for  this  bin.  At  very  early  times  there  is  little 
ionization  due  to  the  ground  source  and  the  200  m  HoB  ionization  rate  is 
much  higher  out  to  about  8  /isec.  Beyond  8  ftsec  the  differences  between 
the  two  sources  are  much  smaller.  In  Fig.  2. 29  a  comparison  between 
the  HoB  =  500  m  and  HoB  =  200  m  sources  for  the  90  -  120  cm  depth 
bin  is  made.  In  this  case,  the  500  m  source  gives  larger  ionization  rates 
at  early  times  to  approximately  5  nsec.  Beyond  5  ^sec  the  ionization 
rate  differences  between  the  two  sources  are  not  very  large  with  the 
HoB  =  200  m  source  ionization  rates  generally  higher. 


Ionization  Rate  (MeV/m  -sec) 


Ionization  Rate  (MeV/m  -sec) 


Figure  2.29.  Ionization  rate  in  ground,  depth  =  90  -  120  cm, 
radius  =  450  -  550  meters,  due  to  thermonuclear 
source  at  HoB  -  200  and  HoB  =  500  meters. 
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A  geometrical  explanation  for  the  differences  in  ground  ioniz¬ 
ation  rates  due  to  source  altitude  variations  is  shown  in  Fig.  2.  30.  The 
point  0  represents  an  observation  point  below  the  air-ground  interface  at 
a  range  of  approximately  500  m.  At  very  early  times  the  ionization  at  0 
is  primarily  due  to  photons  which  are  produced  close  to  the  respective 
neutron  source  locations  through  inelastic  scattering  collisions.  From 
the  figure  it  is  apparent  that  the  prompt  gamma  rays  born  near  the 
HoB  =  500  m  source  must  traverse  less  ground  before  reaching  0  than 
the  prompt  photons  from  the  200  m  HoB.  In  turn,  the  HoB  =  200  m 
prompt  photons  have  far  less  ground  to  traverse  than  gamma  rays  born 
near  the  HoB  =  0  m  source.  Therefore,  at  very  early  times  we  would 
expect  the  ionization  rates  for  the  HoB  =  500  m  source  to  be  larger  than 
the  HoB  =  200  m  rates.  In  addition,  the  HoB  =  200  m  ionization  rates 
at  early  times  would  be  much  larger  than  the  HoB  =  0  m  results  in 
which  the  prompt  photons  are  attenuated  to  a  great  extent  by  the  ground. 

At  later  times,  source  neutrons  have  migrated  away  from  the  respective 
source  locations,  so  that  photons  produced  by  these  neutrons  will  be 
distributed  over  wider  areas.  The  ground  distances  that  must  be  traversed 
by  these  photons  to  reach  point  0,  will  not  differ  to  as  great  an  extent 
between  the  three  sources  as  the  prompt  photons  discussed  above. 
Therefore,  we  would  anticipate  less  deviation  between  ionization  rates 
for  the  three  sources  for  times  later  than  several  microseconds. 

These  remarks  apply  primarily  to  an  observation  point  0  which 
is  deeply  imbedded  in  the  ground,  e.  g. ,  the  bin  discussed  above  at  a 
depth  of  90  - 120  cm.  In  this  case,  the  main  factor  in  determining  the 
ionization  rate  is  the  ground  distance  which  must  be  traversed.  For  bins 
which  are  close  to  the  air-ground  interface,  such  as  the  0-2  cm  depth 
bin  considered  earlier,  photons  can  contribute  to  the  ionization  rate 
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Figure  2.30.  Geometrical  explanation  for  the  effect  of  source 
altitude  on  ground  ionization  rates. 
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without  passing  through  the  ground  to  any  great  degree.  Furthermore, 
direct  neutron  contributions  to  the  ionization  rate  become  important  for 
the?je  bins,  at  local  times  of  several  microseconds,  whereas  for  the 
deeper  bins  the  ionization  is  dominated  by  photons  due  to  their  longer 
mean  free  path  lengths.  Therefore,  the  ionization  rates  for  bins  close 
to  the  air-ground  interface  are  determined  to  some  extent  by  the  amount 
of  air  attenuation  experienced  by  the  photons  and  neutrons  from  the 
various  sources. 

From  these  and  previous  comments,  it  is  clear  that  ionization 
rates  as  well  as  Compton  currents  are  sensitive  to  source  altitudes  in 
some  time  and  spatial  bins  and  relatively  insensitive  in  other  cases.  The 
same  sensitivity  and  insensitivity  can  be  seen  for  other  radial  positions 
as  well.  Arriving  at  general  models  to  interpolate  ionization  rates  and 
currents  between  HoB’s  appears  to  be  a  difficult  task.  Obtaining  analytical 
fits  to  the  data  obtained  for  the  three  HoB's  would  be  a  useful  initial  step 
from  which  possible  interpolation  between  the  various  altitudes  might 
result.  The  effect  of  the  source  altitude  is  a  very  complicated  function 
of  the  spatial  bin  of  interest  as  well  as  the  time  interval  under  consider¬ 
ation.  Further  consideration  of  interpolation  schemes  could  possibly 
result  in  general  models  which  would  eliminate  lengthy  computer  calcu¬ 
lations  for  intermediate  burst  heights. 

In  Figs.  2.  31  and  2.  32  we  also  include  a  brief  comparison  of 
the  present  05RNIES  calculations  and  the  thermonuclear  results  obtained 
last  year.  A  200  m  HoB  is  used  in  both  calculations  with  a  radial  bin 
extending  from  550  -  650  meters  and  an  elevation  bin  extending  from  the 
ground  to  an  altitude  of  50  meters.  Figure  2.  31  compares  ionization  rates 
for  this  bin  while  Fig.  2.  32  shows  the  vertical  current  densities  for  the 
two  calculations.  It  is  difficult  to  make  a  clean  comparison  between  the 

50 


i »-  nki  rtf 


Vertical  Current  Density  (amperes /m 


previous  calculations  and  the  new  results  since  a  new  set  of  neutron  and 
gamma  ray  production  cross  sections  has  been  utilized.  The  sensitivity 
of  transport  and  EMP  calculations  to  cross  section  changes  is  discussed 
more  completely  in  Section  4.  5.  However,  from  Figs.  2.  31  and  2.  32  it 
is  apparent  that  the  inclusion  of  air  capture  events  in  the  new  calculations 
results  in  much  larger  sources  at  very  late  times.  At  early  times,  the 
new  and  previous  results  are  in  reasonable  agreement,  and  we  have  chosen 
to  only  show  the  previous  results  starting  at  approximately  100  /xsec. 
Beyond  100  //sec  the  two  sets  of  calculations  begin  to  deviate  substantially. 
These  large  differences  at  late  times  cannot  be  explained  on  the  basis 
of  cross  section  differences  between  the  two  calculations,  but  must  be 
attributed  to  the  treatment  of  air  capture  reactions  in  the  new  version  of 
05RN1ES.  The  two  comparisons  presented  here  are  representative 
of  calculational  trends  in  other  spatial  bins,  namely  reasonable  agree¬ 
ment  between  the  new  and  previous  calculations  at  early  times,  and  much 
higher  EMP  sources  at  late  times  in  the  new  calculations. 

The  results  presented  to  this  point  have  been  concerned 

primarily  with  the  radial  interval  extending  from  450-  550  meters.  We 

* 

will  now  present  a  number  of  other  radial  intervals  for  the  three  HoB 
sources.  The  results  at  other  ranges  qualitatively  display  the  same 
behavior  discussed  above  for  the  450  -  550  meter  range.  In  Figs.  2.  33  - 
2. 82  we  present  a  number  of  new  calculations  at  ranges  of  850-  1050  m, 
1900-  2100  m,  and  2700-  3100  m  for  the  three  neutron  sources.  We  also 
display  a  few  calculations  for  the  bin  at  4000-  5000  m.  For  this  last 
radial  bin,  the  statistical  uncertainties  in  general  are  large,  since  long 
computer  execution  times  or  extreme  biasing  schemes  are  required  to 
obtain  results  with  small  error  bars. 
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Figure  2.33.  Ionization  rate  at  radius  =  850-  1050  meters, 

cos  6  =  0. 0  -  0. 05,  due  to  thermonuclear  source 
on  ground. 
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Figure  2. 34.  Radial  current  density  vs.  time  at  radius  850 
1050  m,  cos  6  =  0.  0  -  0.  05,  due  to  thermo¬ 
nuclear  source  on  ground. 
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Figure  2. 35.  Theta  current  density  vs.  time  at  radius  =  850  - 
1050  m,  cos  6  =  0. 0  -  0. 05,  due  to  thermonuclear 
source  on  ground. 
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Figure  2.36.  Ionization  rate  vs.  time  at  radius  =  1900  -  2100  m, 
cos  6  =  0. 0  -  0. 05,  due  to  thermonuclear  source 
on  ground. 
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Figure  2.40.  Radial  current  density  vs.  time  at  radius 
1900  -  2100  m,  cos  6  =  0.  5  - 1 . 0,  due  to 
thermonuclear  source  on  ground. 
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Figure  2.46.  Ionization  rate  vs.  time  at  radius  =  4000  - 
5000  m,  cos  0  =  0. 0  -  0. 05,  due  to  thermo 
nuclear  source  on  ground. 
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Figure  2.48.  Ionization  rate  vs.  time  at  radius  =  4000  - 
5000  m,  cos  6  =  0.  5  -  1. 0,  due  to  thermo¬ 
nuclear  source  on  ground. 
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Figure  2.  50.  Ionization  rate  vs.  time  at  radius  =  850  - 

1050  m,  elevation  =  0  -  50  m,  due  to  thermo 
nuclear  source  at  HoB  =  200  m. 
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Figure  2. 51.  Vertical  current  density  vs.  time  at  radius  = 

850  - 1050  m,  elevation  =  0-50m,  due  to  thermo 
nuclear  source  at  HoB  =  200  m. 
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Figure  2.52.  Horizontal  current  density  vs.  time  at  radius 
850-  1050  m,  elevation  =  0-50  m,  due  to 
thermonuclear  source  at  HoB  =  200  m. 
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Figure  2.  55.  Ionization  rate  vs.  time  at  radius  =  1900  - 

2100  m,  elevation  =  0  -  50  m,  due  to  thermo¬ 
nuclear  source  at  HoB  =  200  m. 


76 


Local  Time  (microseconds) 


Vertical  current  density  vs.  time  at  radius  = 
1900  -  2100  m,  elevation  =  0  -  50  m,  due  to 
thermonuclear  source  at  HoB  =  200  m. 
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Figure  2.  60.  Ionization  rate  vs.  time  at  radius  =  2700  - 

3100  m,  elevation  0-50  m,  due  to  thermonuclear 
source  at  HoB  =  200  m. 


Radial  Current  Density  (amperes/na  ) 


Ionization  Rate  (MeV/m  -sec) 


Ionization  Rate  (MeV/m  -sec)  ' 


l  > 


10"2  10’1 


Local  Time  (microseconds) 


Figure  2.  66.  Ionization  rate  vs.  time  at  radius  =  85 
elevation  =  0  -  50  m,  due  to  thermonuc 
at  HoB  =  500  m. 
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Figure  2.  69.  Ionization  rate  vs.  time  at  radius  =  850-1050  m, 
cos  9  =  0.  5  - 1. 0,  due  to  thermonuclear  source  at 
HoB  =  500  m. 
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Figure  2. 70.  Radial  current  density  vs. 

850  -  1050  m,  cos  Q  =  0. 5  - 
nuclear  source  at  HoB  =  5C 
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Figure  2.  76.  Ionization  rate  vs.  time  at  radius  =  2700  -  3100  m 
elevation  0  -  50  m,  due  to  thermonuclear  source 
at  HoB  =  500  m. 
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2700-3100  m,  elevation  =  0-50  m,  due  to 
thermonuclear  source  at  HoB  =  500  m. 
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Figure  2. 78.  Horizontal  current  density  vs.  time  at  radius  = 

2700  -  3100  m,  elevation  =  0  -  50  m,  due  to  thermo¬ 
nuclear  source  at  HoB  =  500  m. 
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Figure  2. 80.  Radial  current  density  vs.  time  at  radius  = 

2700  -3100  m,  cos  =0.5 -1.0,  due  to  thermo¬ 
nuclear  source  at  HoB  =  500  m. 
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Figure  2. 81.  Ionization  rate  vs.  time  at  radius  =  4000- 
5000  m,  cos  6  =  0.  5  -  1.0,  due  to  thermo¬ 
nuclear  source  at  HoB  =  500  m. 
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Ionizaiion  rate  vs.  time  at  radius  =  4000- 
5000  m,  elevation  =  0  -  50  m,  due  to  thermo¬ 
nuclear  source  at  HoB  =  500  m. 
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CURVE  FIT  REPRESENTATION  OF  TRANSPORT  RESULTS 


3. 1  GENERAL  REMARKS 

The  transport  calculations  for  EMP  sources  result  in  a  very 
large  set  of  numbers.  In  addition,  EMP  environment  prediction 
codes  require  sources  at  various  ranges,  angles  and  times  many  of 
which  are  not  directly  available  in  the  transport  calculation  data. 
Therefore,  it  is  necessary  to  provide  some  sort  of  curve  fit  repre¬ 
sentation  to  reduce  the  quantity  of  data  in  the  transport  results  and 
allow  sources  to  be  generated  at  arbitrary  times,  ranges  and  angles. 

This  section  describes  the  curve  fitting  of  transport  data  from 

(3) 

a  previous  EMP  neutron -induced  source  calculation'  '.  This  calcu- 

(4) 

lation  was  for  the  typical  thermonuclear  spectrum'  ,  but  used  an 
older  set  of  the  ENDF/B  data  and  did  not  include  any  air  capture 
effects.  The  transport  bin  structure  is  presented  in  Table  5.  While 
this  previous  calculation  lacks  many  of  the  desirable  features  of  the 
calculation  described  elsewhere  in  this  report,  it  still  makes  available 
a  good  two  dimensional  representation  of  neutron  induced  EMP 
sources. 

The  curve  fitting  effort  was  mainly  directed  toward  the  im¬ 
plementation  of  new  sources  in  the  AFWL  code  SC.  The  new  sources 
will  enhance  the  reliability  of  the  code  in  extended  regions  of  space 
and  time  as  well  as  adding  ground  capture  effects,  angular  dependence 
of  sources  and  transverse  currents  which  were  not  previously  con¬ 
sidered.  While  the  fits  were  produced  with  SC  in  mind  they  should  be 
readily  usable  with  any  other  EMP  ground  burst  code.  They  are 
already  a  part  of  the  AFWL  code  NEW1  and  work  is  proceeding  to 
implement  them  in  SC. 


\ 
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3.2  THE  TRANSPORT  DATA 


The  fitting  effort  considered  only  the  transport  data  in  the 
air  for  the  ground  burst.  Ground  sources  have  been  neglected  as  has 
the  data  for  the  200  meter  height  of  burst.  Aside  from  time  consider¬ 
ations,  the  200  meter  burst  data  was  not  fit  because  there  is  at  present 
no  available  two  dimensional  code  which  could  use  the  data.  One 
dimensional  codes  such  as  NEW1  at  AFWL  can  use  the  data  but  do 
not  seem  to  require  an  elaborate  set  of  fits  in  angle,  range  and  time. 

If  a  two  dimensional  code  is  developed  for  above  ground  bursts,  it 
would  be  desirable  to  base  its  sources  on  the  improved  source  calcu¬ 
lation  described  in  the  other  sections  of  this  report.  Some  studies  with 

NEW1  have  indicated  that  the  ground  sources  associated  with  a  ground 

f  5) 

burst  do  not  change  EMP  fields  drastically  for  close  in  problems. v  ' 
Their  effect  should  be  small  for  ranges  farther  out  and  outside  the 
source  region.  Thus,  since  it  seemed  that  ground  sources  should  not 
be  installed  in  SC  they  were  not  fit.  It  does  seem,  however,  that  the 
fitting  procedure  discussed  here  could  easily  be  made  to  handle  the 
data  for  the  ground  sources  or  the  200  meter  burst  with  minor  modific¬ 
ations  . 

Data  are  available  out  to  a  range  of  3000  meters  and  a  time 
of  10  milliseconds.  However,  statistics  associated  with  the  transport 
results  become  steadily  worse  as  one  goes  to  far  ranges  or  later  times. 
The  fitting  only  considered  times  out  to  1  millisecond.  This  is  appro¬ 
priate  because  air  capture  effects  (  which  were  not  considered)  are 
important  beyond  this  time.  The  statistics  associated  with  the  trans¬ 
verse  (theta)  currents  are  particularly  bad  at  the  far  ranges.  Thus, 
fits  to  the  transverse  currents  only  go  out  to  1000  meters. 
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3.3  THE  FITTING  PROCEDURE 

The  desirable  result  of  a  fitting  effort  is  a  short  FORTRAN 
subroutine  which  can  be  included  in  the  EMP  prediction  codes.  The 
subroutine  must  provide  the  neutron  induced  source  information  for 
arbitrary  points  in  range,  angle  anc  time  and  should  be  reasonably 
efficient  and  easy  to  use. 

It  was  decided  that  it  would  be  exceedingly  difficult  to  find 
functional  forms  which  would  adequately  fit  the  data  in  the  three  vari¬ 
ables  of  interest.  Therefore,  with  the  source  subroutine  in  mind, 
the  following  approach  was  followed.  First,  the  data  were  fit  as  well 
as  possible  in  the  time  domain  for  each  point  in  range  and  angle. 

Then  with  these  time  domain  fits  as  data  for  the  subroutine,  an  inter¬ 
polation  scheme  is  used  to  generate  values  at  any  point  in  range  and 
angle.  Each  time  domain  fit  involves  10  parameters  and  there  are 
285  curves  (5  angle  bins  x  19  range  bins  x  3  source  quantities  =  285 
curves  to  be  fit).  Thus,  the  subroutine  stores  2850  parameters  but 
still  can  occupy  less  than  6000  octal  memory  locations.  This  is  easily 
compatible  with  the  space  available  in  the  EMP  prediction  codes  at 
AFWL.  It  should  be  noted  that  the  first  range  bin  is  deleted  because 
of  poor  statistics  and  that  parameters  are  carried  for  transverse 
current  beyond  1000  meters  though  it  is  recommended  they  not  be  used. 

3.4  THE  TIME  DOMAIN  FITS 

In  the  time  domain  each  source  curve  is  divided  into  three 

regions:  (1)  times  before  neutron  arrival  (based  on  a  15  MeV  neutron); 

(2)  late  times  arbitrarily  defined  to  the  times  greater  than  3  times 

the  neutron  arrival  time;  and  (3)  the  remaining  region  from  arrival 

time  t  to  time  3  t  . 
a  a 
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In  the  region  of  time  tct  the  logarithms  of  the  source  values 
are  fit  as  a  linear  function  of  time  by  a  least  square  polynomial  fitting 
routine  (for  the  transverse  current  the  logarithm  of  1  plus  the  value 
is  fit).  In  the  late-time  region  t  >3t&  a  cubic  function  is  fit  to  the 
logarithm  of  the  source  values  as  a  logarithmic  function  of  time  by  the 
same  routine.  In  the  remaining  region  t  <  t<3t  a  similar  cubic 

ci  cl 

function  is  used  which  provides  continuity  between  the  other  two  fits 
and  satisfies  other  conditions  which  vary  with  the  different  source 

quantities.  Thus,  ten  parameters  Al, - ,  A10  represent  each  curve 

as  follows:  / 

I  exp  (Al  +  A2  t),  t  s  t 

1  a 

Source  (t )  =  (  exp  (A3  +  A4  (In  t)  +  A5  (In  t)*  +  A6  (In  t)3),  t  <  t<  3  t 

1  ci  a 

I  exp  (A7  +  A8  (In  t)  +  A9  (In  t)2  +  A10  (In  t)3),  U3  t  . 

1  ** 

The  parameters  Al,  A2,  A7,  A8,  A9  and  A10  are  calculated 
by  the  least  polynomial  litting  routine  POLY  ^  using  the  transport 
statistics  as  fitting  weight  factors.  The  cubic  function  in  the  middle 
represented  by  parameters  A3,  A4,  A5  and  A6  is  determined  by  the 
following  requirements: 

1.  For  dose  rate  the  function  must 

(a)  connect  with  the  first  fit  at  t  and  the  second  fit 

£L 

at  3t  , 
a  ’ 

(b)  match  the  slope  of  the  second  fit  at  3t  ,  and 

ci 

(c)  pass  through  the  transport  data  point  in  the  first 
time  bin  past 
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2.  For  radial  current  the  function  must 

(a)  connect  with  first  fit  at  t  and  the  second  fit 

cl 

at  3t  ,  and 
a 

(b)  match  the  slopes  of  the  first  fit  at  t  and  the  second 

fit  at  3t  . 
a 

(c)  In  addition,  if  the  resulting  function  crosses  over 

between  t  and  3t  ,  it  is  replaced  by  a  quadratic 

a.  a 

function  which  satisfies  (a)  and  (b)  except  for 

matching  the  slope  at  3t  (in  this  case  A6  =  0). 

cl 

3.  For  transverse  current  the  function  must 

(a)  connect  with  first  fit  at  t  and  second  fit  at 

cl 

3t  and 
a 

(b)  match  the  slopes  of  the  first  fit  at  t  and  the 

cl 

second  fit  at  3t  . 

a 

Examples  of  fits  obtained  in  this  way  are  shown  overlayed 
with  the  transport  data  in  Figs.  3. 1  through  3.26.  It  can  be  seen  that 
the  transverse  current  fits  (JTHETA)  have  some  undesirable  features 
but  the  error  bars  from  the  transport  are  rather  large.  The  iits  beyond 
one  millisecond  should  be  ignored  for  JTHETA.  There  are  some  addi¬ 
tional  difficulties  at  far  ranges  with  the  fits  for  dose  rate  and  radial 
current.  Some  of  them  might  have  been  eliminated  with  additional 
effort  but  it  is  felt  errors  that  are  not  too  drastic  at  far  ranges  will 
not  severely  affect  EMP  field  calculations. 

3.5  THE  INTERPOLATION  PROCEDURE 

To  calculate  source  information  for  use  in  the  EMP  codes, 
interpolation  is  done  between  the  time  domain  fits.  For  dose  rate  and 
radial  current  the  range  interpolation  is  done  linearly  for  the  logarithm 
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of  the  square  of  the  range  times  the  source  quantity.  For  the  trans¬ 
verse  current,  the  range  interpolation  is  linear  for  the  square  of 
the  range  times  the  value.  In  all  cases  the  angle  interpolation  is 
linear  far  the  cosine  of  the  angle. 

For  an  arbitrary  point  in  time,  range  and  angle  (t,  r,  e  ) 
sources  are  calculated  as  follows.  If  the  point  lies  within  the  trans¬ 
port  data  space  grid,  it  will  be  located  between  4  points  where  time 
domain  fits  exist. 


In  the  above  figure  the  four  points  are  labeled  a,  b,  c,  and  d.  Let 
Va,  Vb,  V  Vd  be  the  values  of  any  one  of  the  source  quantities  at 
the  four  points  calculated  from  the  time  domain  fits  at  time  t  .  For 
dose  rate  or  radial  current,  obtain  the  intermediate  quantities  Vj 
and  Vu  by  range  interpolation: 


V,  = 


Vu  ■ 


— z  exP 

r 


~Z  exP 
r 


111  (rl  Vc>  +  F^Fj  (  ln(r2  V  ■  ln  ^1  Vc') 


r  -  r 


(rl  Va>  +  TT71  V  -  to  (rlVa>) 

2  1  2 


For  transverse  current  (which  crosses  over  in  time): 

r  -  r, 


vi  - 


TT 


fr?V+7-Tr;(fr2Vb)-(«-?Vi)) 

Ci  X 


in 


igure  3. 1.  Curve  fits  for  dose  rates  for  typical  thermonuclear 
source  on  ground,  range  =  200  meters,  cos  6  =0.05 


ground. 


Curve  fits  for  theta  current  densities  for  typical  thermo¬ 
nuclear  source  on  ground,  range  =  500  meters,  cos  8  =  0.05 


Local  Time  (microseconds) 

gure  3.7.  Curve  fits  for  dose  rates  for  typical  thermonuclear 
source  on  ground,  range  =  1000  meters,  cos  Q  =  0 


curve  fits  for  radial  current  densities  for  typical  thermo 
nuclear  source  on  ground,  range  =  2000  meters,  cos  6  =  l 


Figure  3.14.  Curve  fits  for  dose  rates  for  typical  thermonuclear 
source  on  ground,  range  =  200  meters,  cos  6  =  0.30. 


nuclear  source  on  ground,  range  =  200  meters,  cos 


.  18.  Curve  fits  for  radial  current  densities  for  typical  thermo¬ 
nuclear  source  on  grouni,  range  =  500  meters,  cos  0  =  0.30. 
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Figure  3.  20.  Curve  fits  for  dose  rates  for  typical  thermonuclear 

source  on  ground,  range  =  1000  meters,  cos  0  =  0.30. 


*T*T» 


curve  fits  for  dose  rates  for  typical  thermonuclear 
source  on  ground,  range  =  2000  meters,  cos  6  =0.30 


2000  MTRS 


Curve  fits  for  radial  current  densities  for  >ical  thermo 
nuclear  source  on  ground,  range  =  2000  meters,  cos  6  =  I 


VU=J7  «r2  V  ‘  (rf  Va» 

r  c  i 


Then  the  value  of  the  source  quantity  is  calculated: 


V  =  V,  + 


COS  Q  -  cos  d. 


1  +  f  COS  $  -  COS  $  1  j  ^u-Vi) 


This  calculation  is  carried  out  in  FORTRAN  routines  which  are 
described  in  Section  3.7. 


3.6  THE  CURVE  FITTING  CODE  NF 

NF  is  a  program  which  reads  a  tape  containing  transport  data, 
performs  the  curve  fits  to  the  data  and  produces  punched  cards  with 
the  curve  fit  parameters  for  use  in  subroutine  NSOURCE.  NSOURCE 
and  its  associated  routines  may  be  used  in  EMP  source  calculations. 

It  performs  the  interpolations  based  on  the  time  domain  fits  and  is 
described  in  the  next  section. 

At  this  time  NF  is  a  somewhat  inefficient  and  poorly  orga¬ 
nized  code.  However,  no  apology  is  made  for  this  situation.  NF  is 
an  assembly  of  many  routines  which  worked  quite  well  in  the  fitting 
effort.  Also,  because  it  has  not  been  made  efficient  or  well  organized 
for  the  effort  described  here,  it  should  be  quite  easy  to  modify  it  for 
use  with  other  data  or  to  perform  other  fitting  tasks. 

A  typical  run  of  NF  takes  only  about  2  minutes  of  computer 
time  (on  the  CDC  6600).  It  requires  the  one  input  tape  with  the  trans¬ 
port  data  and  requires  about  60000  octal  memory  locations. 

A  brief  description  of  NF  follows. 
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3. 7  THE  NEUTRON  INDUCED  SOURCE  SUBROUTINE  NSOURCE 

NSOURCE  has  been  designed  to  calculate  neutron  induced 
source  response  functions  in  a  reasonably  efficient  manner  while  being 
very  easy  to  use.  In  its  present  form  NSOURCE  requires  one  special 
call  which  causes  it  to  read  a  set  of  data  cards.  It  would  be  straight¬ 
forward  for  any  user  to  change  this  situation  to  a  tape  read  if  desired. 

To  use  the  routine  there  must  be  a  one  time  call: 

CALL  NSOURCE  (-1.). 

At  this  time  NSOURCE  will  read  data  cards  and  fill  an  array  of 
arrival  times  for  its  use.  Thereafter  one  may  obtain  the  source 
quantities  with  the  following  call: 

CALL  NSOURCE  (R,  C,  T,  S) 

where 

R  is  the  range  in  meters  where  sources  are  desired, 

C  is  the  cosine  of  the  angle  from  the  vertical, 

T  is  the  time  in  microseconds, 

S  (1)  is  the  returned  dose  rate  value, 

S(2)  is  the  radial  current  value  and 
S(3)  is  the  transverse  current  value. 

A  brief  description  of  NSOURCE  and  its  associated  routines 

follows. 
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Main  Routine: 


NF 


Function:  Organizes  input  data  and  calls  routines  for 

fitting  and  display  of  results. 

Subroutines  called:  AVG,  ARRIVAL,  F,  HISTPLT 

Flow:  1.  Call  AVG  to  fill  arrays  of  transport  bin 

midpoints. 

2.  Read  transport  data  from  TAPE  1. 

3.  Scale  data  if  desired. 

4.  Call  ARRIVAL  to  calculate  neutron  arrival 
times  associated  with  ranges  in  the  transport 
data. 

5.  Call  F  for  fitting  of  each  of  the  three  source 
quantities. 

6.  Call  HISTPLT  to  display  transport  data  and 
fits  on  microfilm  plots. 

7.  Punch  cards  with  fit  parameters. 
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Subroutine  Name: 
Function: 

Called  by: 
Parameters: 


ARRIVAL  (R,  E,  T) 


Calculates  neutron  arrival  time  based 
a  range  and  energy. 

NF 


R  -  range  in  meters 
E  -  neutron  energy  in  MeV 
T  -  neutron  arrival  time  returned 
in  microseconds 


Subroutine  Name: 


F(XX,  YY,  S,  KA,  N,  AB,  AA,  NB,  NA,  NV) 


Function: 

Called  by: 
Subroutines  called: 
Parameters: 


Calls  least  squares  polynomial  fitting  routine 
to  fit  two  parts  of  a  single  source  curve  and 
calculates  and  prints  errors  associated  with 
the  fits. 

NF 

POLY,  POLYF 

XX  -  array  of  times 
YY  -  array  of  source  values 
S  -  array  of  statistics 
KA  -  XX(KA)  closest  point  to  neutron  arrival 
time 

N  -  name  of  quantity  being  fit  e.  g. ,  9  H 
DOSE  RATE 

AB  -  Array  of  parameters  for  fit  before 
arrival 

AA  -  Array  of  parameters  for  fit  after 
arrival 

NB  -  number  of  parameters  in  AB 
NA  -  number  of  parameters  in  AA 
NV  -  variable  flag 
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Subroutine  Name: 
Function: 

Called  by: 


N(X,  F,  A) 

Calculates  parameters  for  a  cubic  function  passing 
through  three  points  (X(l),  F(l)),  (X(2),  F(2)), 
(X(3),  F(3))  with  slope  specified  at  a  fourth  (X(4), 
F(4)). 

PEV 


Parameters: 


X,  F 
A 


see  above 
cubic  parameters 


Subroutine  Name: 


Function: 


POLY  (X,  Y,  Z,  NXYZ,  A,  NA) 

(written  by  Harry  M.  Murphy,  Jr.  of  AFWL) 

Least  squares  polynomial  fit 


Called  by:  F 


arameters:  X  -  array  of  abcissa  values 

Y  -  array  of  ordinate  values 
Z  -  array  of  fitting  weight  values 
NXYZ  -  number  of  values  in  X,  Y  or  Z 

A  -  array  of  polynomial  fit  parameters  returned 
NA  -  number  of  parameters  in  A 
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Function  Name: 


POLYF  (X,  A,  NA) 
POL  YD  (X,  A,  NA) 


Function: 
Called  by: 


Evaluate  a  polynomial  (POLYF)  or  its  derivative 
(POL YD)  at  a  specified  point. 

PEV,  F 


Parameters:  X  -  abcissa  value 

A  -  Array  of  polynomial  coefficients 
NA  -  number  of  values  in  A 
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Subroutine  Name:  PEV(NT,  T,  FT,  TA,  AB,  AA,  NB,  NA,  NE,  QS, 

TS) 

Function:  Evaluate  a  fit  to  a  source  curve 

Called  by:  HISTPLT 

Parameters  NT  -  number  of  times  where  values  are  desired 

T  -  time  array 

FT  -  array  of  fit  values  returned 

TA  -  neutron  arrival  time 

AB  -  array  of  polynomial  fit  parameters  for  curve 
curve  after  arrival 
NB  -  number  of  parameters  in  AB 
NA  -  number  of  parameters  in  AA 
NE  -  source  quantity  flag 

1  -  dose  rate 

2  -  transverse  current 

3  -  radial  current 

(TS,  QS)  point  which  the  dose  rate  fit  is  required  to 
pass  through 


Subroutine  Name: 

Function: 

Called  by: 


AVG(T,  TM,  NT,  R,  RM,  TYPE,  NR  Z 
ZM,  NA,  ND,  NH) 

Fill  arrays  with  transport  bin  midpoint  values 
NF 


Parameters: 


T 

TM 

NT 

R 

RM 

TYPE 

NR 

Z 

ZM 

NA 

ND 

NH 


time  bin  boundaries 

time  bin  midpoints  returned 

number  of  elements  in  TM 

range  bin  boundaries 

range  bin  midpoints 

array  of  flags  for  angle  or  depth 

number  of  elements  in  RM 

cosine  0  or  depth  bin  boundaries 

cosine  0  or  depth  bin  midpoints 

number  of  angles 

number  of  depths 

200  meter  burst  flag 
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Subroutine  Name: 

Function 

Called  by: 
Subroutines  called: 

Parameters: 


,ArZi’  CR’  SCR>  CT>  SCT,  NHORZ 
NVERT,  NTIMES,  RANGE,  CANG  ’ 
HOB,  TNEUT,  TT,  IAN) 


Overlay  fits  on  the  transport  data  for 
plots 


microfilm 


NF 


ix, 


/nn  .  _  7  *^wxvxiJXJX\,  STKjV 

?YMB0L  AND  NUMBER  are  parts  of  a 
standard  plot  package  of  AFWL) 


Q 

SQ 

CR 

SCR 

CT 

SCT 

NHORZ 

NVERT 

NTIMES 

RANGE 

CANG 

HOB 

TNEUT 

TT 

IAN 


transport  dose  rate  values 
statistics  associated  with  dose  rate  values 
transport  radial  current  values 
statistics  associated  with  radial  current 
values 

transport  transverse  current  values 
statistics  associated  with  transverse 
current  values 

number  of  decades  horizontally  in  plot 
number  of  decades  vertically  in  plot 
number  of  elements  in  Q,  SQ  etc 
range 

value  of  cosine  0 
height  of  burst 
neutron  arrival  time 
transport  time  bin  mid  points 
number  of  the  angle  bin 
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Subroutine  Name:  QUAD  (Yl,  Tl,  Y2,  T2,  YPZ,  A,  B.  C) 


Function: 


Called  by: 


Calculates  parameters  for  a  quadratic  function 
passing  through  two  points  with  slope  specified 
at  the  third. 

PEV 


Parameters 


(Tl,  Yl) 

(T2  Y2)  "  Two  sPecifie<*  points 

YP2  -  specified  slope  at  T2 
A,  B,  C  -  quadratic  parameters 
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Subroutine  Name: 

Function: 

Called  by: 
Parameters: 


'i"iy  ' 


'«w^l 


SPLINE  (Yl,  YP1,  Tl,  Y2,  YP2,  T2,  A,  B.  C  D 
IDRIV) 


Calculates  parameters  for  a  cubic  function  passing 
through  two  points  with  slope  specified  at  each. 

PEV 


Two  specified  points 


(Tl,  Yl) 

(T2,  Ye) 

YP1,  YP2  slopes  at  Tl  and  T2 

A,  B,  C,  D-cubic  parameters 


IDRIV  -flag  which  indicates  whether  or  not 

calculated  cubic  crosses  over  between 
Tl  and  T2. 


150 


Subroutine  Name:  NSOURCE  (RX,  CX.  TX,  SX) 

Function:  To  calculate  neutron  induced  EMP  source  terms  at 

any  range  angle  and  time  covered  by  the  transport 
data. 

Subprograms  called:  NSET,  VAL,  XINT,  XINTL 

Flow:  1.  Check  to  see  if  RX  is  negative. 

If  so,  call  NSET  and  return. 

2.  Locate  (RX,  CX)  in  mesh  of  points  where 
time  domain  fits  are  available 

3.  Call  VAL,  XINT  and  XINTL  to  interpolate 
for  source  term  values. 


Function  Name: 


Function: 

Called  by: 
Parameters 


XINTL  (R,  C,  RA,  RB,  CA,  CB,  VA  VB  VC  VDl 
XINT  (R,  Cf  RA,  RB,  CA,  CB,  VA.  VBVCW) 


Perform  the  linear  interpolation  for  transverse 
current  value  (XINTL)  of  the  log  interpolation  for 
dose  rate  or  radial  current  (XINT) 


NSOURCE 


R  -  r  (See  Section  3,4) 
C  —  c 
RA  -  r, 

RB  -  ri 
CA  -  cf 
CB  -  co 
VA  -  v 
VB  -  V* 

VC  -  Vb 

VD  -  vf 
d 
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Subroutine  Name: 
Function: 

Called  by: 
Parameters 


NSET  (NR,  R,  TA,  DD.  DR,  DT) 

To  read  time  domain  fit  parameters  and  fill 
arrival  time  array  for  later  use  in  NSOURCE 

NSOURCE 

NR  -  number  of  ranges  where  time  domain  fits 
are  available 
R  -  array  of  ranges 
TA  -  array  of  arrival  times 
DD  -  array  of  time  domain  fit  parameters  for 
dose  rate 

DR  -  array  or  time  domain  fit  parameters  for 
radial  current 

DT  -  array  of  time  domain  fit  parameters  for 
transverse  current 


mmrn 


•w  ji  1  j  ^  hw^jw.%  ■■  ^v”T  jfy  v* 


Function  name: 
Function: 

Called  by: 
Parameters: 


VAL  (T,  N,  A) 

to  evaluate  a  time  domain  fit 

NSOURCE 

T  -  time  for  evaluation 

N  -  =  1  use  fit  for  t  <  t 

—  a 

=  2  use  fit  for  t  <  t  <  3t 
a  a 

=  3  use  fit  for  t  >  3t 
a—  a 

A  -  array  of  fit  parameters 
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4.  THE  05RNIES  PROGRAM 


4. 1  GENERAL  REMARKS 

The  05RNIES  program  is  a  versatile  Monte  Carlo  code  for 

solving  the  time-dependent  neutron  and  gamma  ray  transport  equations 

in  two  space  dimensions.  From  these  solutions,  the  Compton  currents 

and  ionization  rates  in  air  and  ground  due  to  a  specified  neutron  source 

are  obtained.  Late-time  electromagnetic  fields  may  then  be  computed 

by  solving  Maxwell’s  equations  together  with  the  relevant  air  and  ground 

chemistry  equations.  A  brief  discussion  of  these  equations  was  presented 

m 

in  a  previous  report  and  will  not  be  repeated  here. '  '  The  05RNIES 
program  described  in  this  section  has  been  substantially  revised  from 
the  version  of  last  year.  We  will  attempt  to  emphasize  those  areas  of 
the  code  which  have  been  altered  to  any  significant  extent.  We  concentrate 
here  on  the  newer  portions  of  the  code  and  include  some  of  the  revised 
subroutines  of  last  year. 

The  new  version  of  05RNIES  has  been  programmed  with  a  band 
source  energy  concept  to  facilitate  folding  in  arbitrary  neutron  spectra 
at  a  later  time.  In  addition,  new  air  capture  techniques  have  been 
introduced  which  extend  the  time  scale  out  to  the  one  second  region. 
Neutrons  in  the  lower  energy  ranges  are  now  tracked  from  one  medium 
to  another  so  that  capture  gamma  rays  in  both  air  and  ground  are 
accurately  taken  into  account.  Thermal  neutron  energy  deposition  in  the 
air  is  also  considered  for  the  first  time.  New  gamma  ray  analytical 
methods  have  also  been  introduced  for  transporting  gamma  rays  at 
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positions  distant  from  the  air -ground  interface.  Individual  contributions 
to  ionization  and  Compton  currents  by  neutrons  and  gamma  rays,  as  well 
as  by  fast  and  thermal  reaction  processes,  are  now  displayed  to  assist  in 
fitting  of  the  vast  amount  of  data  generated. 

The  latest  ENDF/B  neutron  transport  and  gamma  ray  production 
cross  sections  have  been  processed  and  incorporated  into  data  files  for 
use  in  these  calculations.  Ionization  rates  and  Compton  currents  are 
stored  on  magnetic  tapes  at  AFWL  as  these  data  are  produced  so  that 
they  may  be  retrieved  for  later  calculations.  05RNIES  also  has  restart 
capabilities  which  enable  the  user  to  finish  a  calculation  which  was 
aborted  for  any  reason. 

4. 2  SOURCE  ENERGY  BAND  VERSION 

4.2.1  Selection  of  Bands 

Full  production  calculations  using  05RNIES  may  typically  require 
20  hours  of  time  on  the  CDC  6600  to  generate  200,000  neutron  histories. 

In  the  past,  a  new  05RNIES  calculation  was  needed  for  each  input  neutron 
source  spectrum.  Since  there  are  many  neutron  spectra  of  interest  in 
EMP  studies,  it  became  clear  that  a  new  computational  scheme,  which 
would  not  require  rerunning  full  production  calculations  for  each  new 
neutron  source,  was  needed.  The  source  energy  band  version  of  05RNIES 
was  designed  to  fulfill  this  requirement. 

The  present  version  of  the  code  has  a  total  of  13  fixed  energy 
bands.  The  number  13  was  chosen  as  a  compromise  between  the  many 
groups  actually  required  to  accurately  fit  arbitrary  nuclear  weapon  neutron 
output  spectra  and  only  a  few  groups  which  would  reduce  the  computer 
execution  times.  The  new  code  is  run  until  a  sufficient  number  of 
neutron  histories  are  generated  for  each  source  band.  Various 
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quantities  of  interest  are  then  stored  on  data  tapes  as  a  function  of  neu¬ 
tron  source  group.  The  auxiliary  folding  code  05NPT  is  then  used  to 
generate  EMP  drivers  for  any  desired  input  neutron  source  spectrum. 
In  this  way,  only  a  few  minutes  of  computer  time  are  required  to  gen¬ 
erate  EMP  sources  for  a  new  neutron  source  rather  than  the  20-hour 
calculations  of  previous  years. 

In  Table  6  the  13  energy  bands  used  in  these  calculations  sre 
displayed.  The  probability  distribution  of  source  neutrons  for  the  typi¬ 
cal  thermonuclear  spectrum  of  Table  1  is  also  shown  in  the  new  13- 
group  structure.  In  Fig.  4. 1  a  comparison  of  the  new  13-band  source 
distribution  and  the  original  17-group  thermonuclear  distribution  is 
made.  It  is  apparent  that  the  new  band  structure  is  very  satisfactory 
for  the  thermonuclear  spectrum.  Corroboration  of  the  adequacy  of  the 
band  structure  for  several  other  devices  has  also  been  made. 

4.2.2  Statistical  Analysis 

Estimates  for  the  mean  and  variance  of  computed  quantities 
using  the  new  band  version  of  05RNIES  are  obtained  in  a  straightfor¬ 
ward  fashion.  Assume  that  there  are  M  energy  bands  of  source  neu¬ 
trons,  and  for  each  band  there  are  n  neutron  histories  followed  in 
each  batch  or  run.  There  are  then  N  batches  or  rims  performed  to 
generate  production  results  for  a  given  HoB.  Furthermore,  let  ffc 
be  the  fraction  of  source  neutrons  in  the  k-th  energy  band.  Then  if 
X...  represents  the  i-th  batch  sample,  for  history  index  j  and  band  k, 

1JK 

we  easily  find  for  the  i-th  sample, 


M  n 
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Table  6 


COMPARISON  OF  THE  TYPICAL  THERMONUCLEAR 
NEUTRON  SOURCE  SPECTRUM  IN  THE 
ORIGINAL  17 -GROUP  STRUCTURE  AND  IN  THE 
NEW  13 -BAND  SOURCE  VERSION. 


17 -Group  Thermonuclear 
Neutron  Distribution 

Band 

No. 

13 -Band  Energy 

Version 

Group 

No. 

Upper  Energy 
(MeV) 

Probability 

Upper  Energy 
(MeV) 

Probability 

1 

15.000000 

. 07060 

1 

16.500000 

. 00000 

2 

12.200000 

.02560 

2 

15.000000 

. 03782 

3 

10.000000 

.01410 

3 

13.  500000 

. 06578 

4 

8.190000 

. 01470 

4 

9.  050000 

.03959 

5 

6.360000 

.01800 

5 

4.  960000 

. 06681 

6 

4.970000 

. 01700 

6 

2.350000 

.09152 

7 

4.070000 

. 02600 

7 

1.100000 

. 08348 

8 

3.010000 

.01900 

8 

. 550000 

.07186 

9 

2.460000 

. 00500 

9 

. 240000 

.  03014 

10 

2.350000 

. 02800 

10 

.110000 

.26832 

11 

1.830000 

. 06200 

11 

.031600 

.04210 

12 

1.110000 

.  08500 

12 

.019300 

.  05455 

13 

. 550000 

.10200 

13 

. 003360 

.  14803 

14 

.110000 

.36500 

. 000023 

15 

. 003350 

.  12200 

16 

. 000583 

.  02400 

17 

. 000101 

.  00200 

.000029 

158 


Probab  il  ity/  Me  V 


Energy  (MeV) 


Figure  4. 1.  Comparison  of  the  typical  thermonuclear  spectrum 
in  the  original  17-group  structure  and  in  the  new 
13-band  version. 
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Then  we  also  have 


The  estimates  for  the  mean  and  variance  are  then  found  from 


and 


"2=  ~!~rr  £  X2 
NtfPIJ  Xi  - 

In  our  sampling  we  perform  N  runs  or  batches,  each  run  con¬ 
sisting  of  M  .  n  neutron  histories.  The  unit  sample  is  the  run  and  not 
an  individual  neutron  history.  The  principal  advantage  to  this  method 
of  group  mg  is  the  reduction  in  the  number  of  samples  and  sample 
squares  which  must  be  computed  and  stored  for  the  statistical  analysis. 
In  our  calculations  at  each  of  three  HoB's,  we  have  performed  a  total 
of  N  =  15  runs,  with  M  =  13  energy  bands,  and  n  =  1000  neutron  his¬ 
tories  for  each  band.  Thus,  each  production  run  is  comprised  of  a 
total  of  195,000  neutron  histories. 


4.  3  LATE- TIME  EMP  DRIVERS 

4.3.1  Air  Capture  Computational  Models 


(1-2) 

In  earlier  versions  of  05RNIES'  '  there  were  limitations  on 
the  validity  of  results  of  late  times  because  contributions  from  thermal 
neutron  capture  in  the  ground  were  not  considered.  Last  year  a  signif¬ 
icant  improvement  in  the  05RNIES  late-time  physics  calculations  was 
made  by  introducing  a  time-dependent  one- velocity  diffusion  model  to 
treat  thermal  neutron  capture  events  in  the  ground.  The  model,  although 

4 

approximate,  enacted  calculations  out  to  times  as  late  as  10  ^sec.  The 
diffusion  model  rei  taced  the  air,  in  the  two-region,  air-ground  problem, 
by  a  vacuum,  so  that  thermal  neutrons  leaving  the  ground  and  entering 
the  air  never  returned  to  the  ground.  This  approximation  is  valid  for 
the  time  range  discussed  above.  However,  when  air  capture  events  are 

a 

considered,  local  times  as  late  as  10  nsec  are  of  interest  and  a  new 
computational  model  which  correctly  tracks  low  energy  neutrons  from 
one  medium  to  another  is  required.  The  air  capture  models  described 
below  are  applicable  to  both  air  and  ground  regions  and  enable  accurate 
calculations  to  times  of  1  second. 


The  two-region,  time-dependent,  air  capture  gamma  problem 
has  been  one  of  the  largely  unresolved  radiation  transport  problems  of 
importance  to  EMP  studies  at  late  times  (beyond  about  10  f/sec).  A  direct 
brute  force  Monte  Carlo  analysis  of  this  problem  was  not  computationally 
efficient  because  of  the  extremely  long  computer  times  required.  The 
difficulty  stems  mainly  from  the  fact  that  air  capture  events  take  place 
predominantly  at  low  neutron  energies  and  many  neutron  collisions  are 
required  before  thermalization  in  air  is  achieved.  For  example,  a 


161 


; 


neutron  in  air  suffers  approximately  five  times  as  many  collisions  in 
downscattering  from  0. 1  MeV  to  thermal  energies  as  it  does  in  slowing 
down  from  a  typical  source  energy  of  2.  0  MeV  to  0. 1  MeV.  Therefore, 
direct  Monte  Carlo  calculations  in  the  epithermal  range  (0. 1  MeV  to 
2.  38  eV)  as  well  as  in  the  thermal  group  (below  2.  38  eV)  could  con¬ 
ceivably  take  several  times  as  long  to  perform  as  the  previous  25 -hour 

(3) 

runs  required  for  fast  neutron  EMP  studies.  Biasing  of  the  Monte  Carlo 
calculations  to  emphasize  air  capture  events  would  be  one  way  to  lower 
computer  execution  times  for  this  particular  problem.  Some  efforts 
along  these  lines  were  pursued,  but  eventually  abandoned  in  favor  of 
new  analytical  approaches  which  show  great  promise. 

In  the  calculations  discussed  here  the  neutron  transport  is 

carried  out  by  different  methods  in  each  of  three  separate  energy 

ranges.  In  the  highest  energy  range,  above  an  input  quantity  E  , 

c 

usually  about  0.  Ill  MeV,  the  neutron  transport  is  analyzed  by  the 
standard  Monte  Carlo  approach.  In  the  epithermal  range  between 
E  and  the  thermal  group  cutoff  energy,  E.  =  2.38  eV,  an  analytical 
age  diffusion  method  is  coupled  with  the  conventional  Monte  Carlo 
technique.  The  analytical  model  is  based  on  the  premise  that  a  neu¬ 
tron  which  is  sufficiently  high  above  the  air- ground  interface  can  be 
accurately  treated  by  infinite  medium  Fermi  age  theory.  Similarly, 
neutrons  which  are  in  the  ground  and  sufficiently  below  the  air -ground 
interface  can  also  be  treated  by  age  theory.  The  analytical  scheme 
is  used  for  spatial  positions  distant  from  the  air-ground  interface, 
whereas  the  Monte  Carlo  approach  is  used  for  particles  near  the  two- 
region  interface. 
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For  thermal  neutrons,  E  <  E^,  a  one-group,  time-dependent, 
diffusion  model  is  coupled  with  a  one-velocity  Monte  Carlo  technique. 

The  diffusion  theory  approach  can  be  employed  whenever  thermal  neu¬ 
trons  are  located  at  spatial  positions  distant  from  the  air-ground  inter¬ 
face,  while  the  Monte  Carlo  method  is  used  near  the  interface.  When 
a  thermal  neutron  capture  gamma  ray  is  produced,  a  new  analytical 
gamma  ray  transport  model  is  also  introduced  to  track  those  photons, 
in  both  air  and  ground,  whic..  are  initially  distant  from  the  air-ground 
interface. 

by  coupling  of  the  Monte  Carlo  and  the  analytical  methods,  the 
approximate  analytical  techniques  are  applied,  where  appropriate  in 
the  epithermal  and  thermal  groups,  to  greatly  improve  computer  exe¬ 
cution  times.  On  the  other  hand,  the  Monte  Carlo  method  is  retained 
to  treat  particles  in  spatial  or  energy  regimes  in  which  the  analytical 
methods  are  not  valid.  The  coupled  techniques  therefore  enable  accur¬ 
ate  extension  of  previous  transport  calculations  to  later  time  regions 
with  good  computational  efficiency. 

In  the  next  three  sections  we  describe  in  more  detail  how  epi¬ 
thermal  and  thermal  neutron  transport,  as  well  as  gamma  ray  transport, 
are  dealt  with  in  the  new  code. 

4.3.2  Epithermal  Neutron  Transport 

Since  a  brute  force  approach  to  solving  the  late-time  air  cap¬ 
ture  problem  was  not  deemed  satisfactory,  an  alternative  analytical 
procedure  was  developed.  The  new  model  basically  determines  whether 
or  not  a  neutron  in  the  air  is  sufficiently  distant  from  the  ground  so 
that  simpler  transport  approximations  can  be  incorporated.  For  neu¬ 
trons  far  from  the  air-ground  interface,  the  new  approximations 
assume  that  complexities  introduced  by  the  ground  can  essentially  be 
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ignored  and  air  transport  calculations  can  be  performed  with  an  infinite 
medium  Fermi  age  model.  Similarly,  for  neutrons  deeply  imbedded  in 
the  ground,  infinite  medium  ground  transport  calculations  are  assumed 
to  be  appropriate  with  the  air  being  ignored.  In  this  section  we  describe 
the  coupled  Monte  Carlo-age  theory  calculations  which  are  used  in  treat¬ 
ing  neutrons  in  the  epithermal  energy  region  E  >  E  >  E, .  When  neu- 

c  t 

trons  are  downscattered  to  energies  below  E^,  one  velocity  Monte  Carlo 
calculations  are  combined  with  time-dependent  thermal  diffusion  theory 
calculations  as  described  in  the  next  section. 


For  a  given  neutron  of  energy  E  and  altitude  Z  above  the  air- 
ground  interface,  it  is  possible  from  age  theory  and  solid  angle  argu¬ 
ments  to  compute  the  probability,  P(Z,  E-E^),  that  the  neutron  in 
slowing  down  to  the  thermal  energy  group  cutoff,  Et,  will  strike  the 
ground.  If  P(Z,  E— Et)  is  smaller  than  some  small  predetermined 
lower  probability  limit,  PM^x>  then  the  two-region,  air-ground  prob¬ 
lem  can  be  effectively  uncoupled,  and  the  neutron  can  be  treated  by 
infinite-air  Fermi  age  theory.  If  P(Z,  E— Et)  is  not  smaller  than  Pj^j^> 
a  new  lower  energy  E^  may  be  selected,  where  E^  >  Et,  so  that 
P(Z,  E— Ej)  is  smaller  than  P^j^-  Age  theory  infinite  medium  cal¬ 
culations  can  then  again  be  applied. 

Eventually,  there  is  an  altitude,  Z  ,  below  which  it  is  not 

c 

appropriate  to  employ  infinite  medium  age  theory  since  there  is  an 

intolerably  high  probability  for  a  neutron  to  cross  from  the  air  to  the 

ground  for  this  altitude.  Note  that  this  is  not  the  case  if  we  chose  to 

only  degrade  the  neutron  an  infinitesimal  amount  in  energy,  E.— >E.  -A 

J  J 

where  A— *0.  However,  this  would  not  be  efficient  computationally 
and  could  also  introduce  a  bias  into  the  calculations.  Therefore,  in 
our  calculations  for  distances  in  air  below  Z  =  75  meters,  we  have 
used  the  Monte  Carlo  method  exclusively. 


The  Monte  Carl°  ‘^hnique  is  used  to  follow  the  neutron  until 

‘  e  ^  e,haS  6ither  lnCreafied  ab0TC  ”  meters  or  until  it  has  reached 
he  thermal  neutron  group,  E  E  E(.  As  the  neutron  moves  away  from 

the  air-ground  interface  it  is  monitored  and  analytical  methods  are  used 

as  soon  as  appropriate.  A  neutron  may  be  effectively  followed  by  an 

alternating  succession  of  Monte  Carlo  and  analytical  routines  if  the 

neu  ron  travels  back  and  forth  between  the  narrow  Monte  Carlo  region 

“r‘0nrWay  'r0m  the  air'gr°“d  in‘erface-  The  comhined  Monte 
Carlo-age  theory  tracking  of  the  neutron  is  continued  until  the  neutron 

energy  falls  within  the  thermal  energy  group 


Irons  '  it”  entlrely  Slmilar  process  is  used  to  analytically  treat  neu¬ 
rons  in  the  epithermal  range  which  are  in  the  ground.  In  this  case  the 

analytical  method  is  used  for  depths  greater  than  5  cm  and  the  Monte 

ar  o  technique  is  used  to  follow  neutrons  which  are  close  to  the  air- 
ground  interface. 


The  neutron  flux  in  the  epithermal  energy  range,  at  spatial 
positions  distant  from  the  air-ground  interface,  is  determined  by  solv- 
g  he  age-diffusion  equation.  For  a  monoenergetic  spherically 
ymmetric  point  source  at  the  origin  the  age  equation  may  be  written  as 


D($  +  f  v  + 

4  nr 


S  (t)  5(u)5  (4.1) 


where 


$  =  0(r,  u,  t) 


(4.2) 
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and  the  lethargy  u  is  related  to  the  neutron  speed  v  by 


v  = 


(4.3) 


where  vQ  is  the  neutron  source  speed.  The  other  notation  in  Eq.(4. 1  ) 
is  standard. 

For  an  infinite  medium,  Eq.  (4. 1)  is  solved  by  applying  succes¬ 
sively  Fourier  spatial  and  Laplace  time  transforms  and  then  inverting. 

The  solution  is  found  to  be 

,  ,  Qo  „-(r2/4  r) 

*<r- »■ ' l>  -  fyr e  J^37r  ^  (t-‘sD(Vv»  >  <4- 4) 

where  the  slowing  down  time  tgD  is  defined  by 


(v-v) 


fv° 

1  i 


s 


(4.5) 


and  the  Fermi  age  r  is  given  by 


r(0— u) 


r(Eo-E) 


D[u')da’  _  [E° D(E’)  dE’ 

J£  fX^E'  " 


(4.6) 
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To  be  able  to  use  the  neutron  probability  distribution  functions  of 
Eq.  (4. 4),  an  expression  must  be  derived  which  relates  the  validity  of 
the  infinite  medium  age  equation  to  the  altitude,  or  depth,  Z,  and  the 
incident  energy  E  of  the  neutron. 

We  concern  ourselves  henceforth  with  the  case  of  a  neutron  in 
the  air;  the  problem  of  a  neutron  in  the  ground  is  treated  by  very  simi¬ 
lar  methods.  We  have  found  numerically  that  the  air  and  ground  can 
be  effectively  uncoupled  and  age  theory  can  be  introduced  if  the  altitude 
of  the  neutron  is  greater  than 


Z  =  C  V(r(H^-Ea) 


(4.7) 


Here  Z  is  the  neutron  altitude  above  the  air-ground  interface  and  C 
is  a  constant.  Values  of  C  between  2.  5  and  3.  5  have  been  found  to 
give  satisfactory  numerical  results  in  our  calculations.  And  finally, 
T  Ea)  is  the  Fermi  age  for  downscattering  from  initial  energy 
E^  to  final  energy  Ea 


% 


■E  )  - 
a 


■/' 

JEr 


*b  D(E') 


OglT) 


dE’ 

E’ 


(4.8) 


An  expression  equivalent  to  Eq.  (4.7)  holds  for  a  neutron  in  the  ground 
where  Z  in  that  case  refers  to  the  depth  below  the  interface. 


Physically  r(E --E  )  of  Eq.  (4.  8)  is  a  measure  of  the  mean  square 

D  3 

distance  traveled  by  a  neutron  as  it  downscatters  from  energy  E^  to 
energy  Ea.  Hence,  as  Ea  is  reduced,  r(Eg-Ea)  and  the  distance 
traveled  by  the  neutron  increases,  and  the  altitude  Z  of  Eq.  (4.7)  must 
also  increase  if  the  infinite  medium  theory  is  to  be  valid. 
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We  are  given  an  epithermal  neutron  of  energy  and  altitude 
Z  above  the  ground.  The  problem  is  to  determine  the  final  energy  E 

I cl 

in  order  that  age  theory  hold.  It  is  apparent  that  as  E  approaches  E.  , 
r(EjpEa)  of  Eq.  (4.  8)  tends  to  zero  and  a  very  low  altitude  Z-*0  will 
satisfy  Eq.  (4.7).  Values  of  E  very  close  to  E,  are  not  satisfactory 
from  a  computational  standpoint  since  bias  can  be  introduced  into  the 
results.  Therefore,  we  only  allow  age  theory  calculations  to  take  place 
if  a  non-negligible  energy  decrement,  AE,  of  the  order  of  10%  is  avail¬ 
able.  We  furthermore  only  permit  age  theory  calculations  in  the  air 
above  75  meters  and  in  the  ground  at  depths  below  5  centimeters.  In 
this  way  Monte  Carlo  is  used  to  treat  the  important  region  near  the  air- 
ground  interface  without  introducing  a  bias  due  to  overextension  of  the 
age  theory  models. 

We  use  Eqs.  (4. 7)  and  (4.  8)  in  the  following  manner  for  an  epi¬ 
thermal  neutron  in  air.  An  initial  test  is  made  to  see  if  the  neutron  is  above 
the  75-meter  Monte  Carlo  cutoff  plane.  The  final  energy  E  of  the  neutron, 

ci 

after  downscattering  through  a  number  of  scattering  collisions,  is  then 
determined  by  an  iterative  scheme.  We  would  like  E  to  be  as  small  an 
energy  as  possible,  preferably  E^.,  so  that  the  normally  time  consuming 
Monte  Carlo  downscattering  process  in  the  epithermal  range  is  replaced 
by  the  efficient  analytical  age  theory  models.  It  is  clear  that  if  the  neutron 
is  very  high  above  the  ground,  it  should  be  possible  to  downscatter  the 
neutron  to  the  thermal  group  cutoff,  E  =  E.,  analytically.  For  example, 
if  a  neutron  of  energy  0.  Ill  MeV  is  located  700  meters  above  the 
ground,  the  probability  of  analytically  thermalizing  the  particle  is 
very  high.  In  this  case,  the  probability  of  the  neutron  reaching  Et  by 
age  theory  calculations,  without  striking  the  ground,  is  greater  than 
0.  99.  On  the  other  hand,  a  0.  Ill  MeV  neutron  at  an  altitude  of  100 
meters  has  a  much  higher  probability  of  striking  the  ground  if  it  is 
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allowed  to  downscatter  to  thermal  energies.  Therefore,  for  the  lower 
altitude  particle,  smaller  energy  decrements  are  allowed  and  new  spatial 
positions  are  computed  by  age  theory.  If  the  new  altitudes  which  are 
computed  for  the  particle  increase  over  the  initial  100-meter  height,  the 
probability  of  eventually  thermalizing  the  neutron  by  age  theory  calcula¬ 
tions  is  also  increased.  However,  if  the  particle  drops  below  the  75- 
meter  cutoff  plane,  the  conventional  Monte  Carlo  method  is  again  insti¬ 
tuted  until  the  particle  either  is  scattered  above  75  meters,  or  until  it 
is  downscattered  and  reaches  thermal  energies.  Another  possibility 
is  for  a  Monte  Carlo  neutron  in  the  air  to  be  tracked  deep  into  the  ground 
and  eventually  to  reach  such  a  depth  that  a  ground  age  theory  calcula¬ 
tion  may  be  used. 

The  calculation  of  the  final  allowable  energy  E  is  derived 

3. 

from  squaring  of  Eq.  (4.7), 


Z 


2 


„2  ...  fE b  D(E')  dE'  _  „2 

JE  njm  *•  - 


(VEa> 


(4.9) 


or  on  inverting  Eq.  (4.9), 


Ea  =  Ejj  x  exp 


(Z/C)2 

<D/f£s) 


(4. 10) 


where  (D/^£  )  is  a  weakly  varying  function  of  energy  and  is  averaged 
over  the  energy  range  between  E,  and  E  .  From  Eq.  (4. 10)  it  is 
apparent  that  once  the  altitude  Z  and  incident  energy  are  known, 
it  is  possible  to  iteratively  determine  the  final  energy  E  after  energy 

a 

averaging  the  (D/^2  )  cross  section  terms  which  depend  on  E,  and  E  . 

S  D  cl 

Once  E  is  found,  the  age  r(E,  E  )  can  also  be  computed  from 
Eqs.  (4.  8)  and  (4.  £).  This  age  then  represents  the  largest  age  decrement 
which  a  neutron,  at  altitude  Z  with  initial  energy  Eb,  can  undergo  with¬ 
out  accepting  an  intolerably  high  probability  of  striking  the  ground. 

The  actual  spatial  position  of  the  neutron  following  this  aging  process 
is  found  by  sampling  a  distance  r  from  a  time- integrated  and  normal- 
lized  solution  of  the  age  equation  of  Eq.  (4.  4).  This  solution  or  prob¬ 
ability  density  function  is  given  by 

e"1*2  /4r  .  2 

i//(r,  u)  =  t/f(r,  r )  = - ^  4ur,  (4.11) 

(4  »7  r  )  i/l 

where  the  r  in  Eq.  (4.  11)  is  the  r(Ec-*E  )  as  determined  by  Eqs.  (4.  8) 

D  a  , 

and  (4.  9).  The  distance  r  traveled  by  a  neutron  in  aging  from  E^  to 
Ea  is  found  by  sampling  from  Eq.  (4.  11)  by  a  rejection  technique.  The 
initial  spatial  coordinates  (Xj,  y^,  z^)  as  well  as  the  initial  time  tj 
are  known.  To  find  the  final  space  coordinates,  we  generate  an  isotropic 
unit  vector  fl  with  components  a,  / 3,  and  yin  the  x,  y,  and  z  directions. 

The  final  spatial  position  is  then  (Xg,  where 

Xg  =  Xj  +  r  .  a 

y2  =  yx  +  r  .  /3  (4.  12) 

z2  =  z1+  r  •  y 

where  r  is  sampled  from  Eq.  (4. 11). 
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In  the  relatively  rare  case  where  the  isotropic  vector  Q  points 
straight  downwards,  it  is  possible  for  a  neutron  in  the  air  to  cross  over 
into  the  ground.  In  these  exceptional  cases,  a  new  isotropic  vector 
is  selected  so  that  the  particle  does  not  move  from  one  media  to  another. 
These  selection  processes  have  been  carefully  verified  so  that  essentially 
no  bias  is  introduced  into  the  calculations. 

The  time  required  for  the  neutron  to  age  from  energy  to 

energy  E  is  given  by  the  slowing  down  time  t__  (v,-*v  )  of  Eq.  (4.  5). 

The  total  time,  from  birth  as  a  source  neutron  to  energy  E  is  given  by 

21 


42  "  ll  +  tSD(VVa), 


(4. 13) 


where  t^  is  the  net  time  to  the  particles  appearance  at  energy  E^  and 
tgD  is  the  slowing  down  time  to  energy  Ea* 

We  then  have  the  final  particle  position  (x2,  y2,  z2)  and  time  t2 
which  become  the  initial  coordinates  for  the  next  downscattering  or 


thermal  diffusion  calculations. 


In  Eq.  (4. 11)  it  is  apparent  that  neutron  absorption  is  not  properly 
considered.  The  probability  of  neutron  absorption  during  the  aging  process, 
E^pE  ,  is  included  by  using  the  concept  of  a  resonance  escape  probability 
p(EtrEa),  where 


P(Eb-rEa>  =  e*P  j-2a(v)vWVEa)} 


(4. 14) 


The  neutron  weight  before  collision,  W^,  is  multiplied  by  p(E^-*Ea)  to 


obtain  the  weight  after  aging,  W  , 

cl 


w„  =  w„  p(E.— E„ ). 


(4. 15) 
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This  is  analagous  to  the  usual  Monte  Carlo  treatment  of  absorp¬ 
tion  where  the  neutron  weight  is  multiplied  by  the  nonabsorption  proba¬ 
bility,  (Sg/Ij.),  to  determine  the  final  weight. 

To  determine  a  spatial  position  and  time  for  depositing  energy 
due  to  neutron  absorption  and  elastic  scattering  in  the  epithermal  energy 
range,  we  first  generate  a  new  age  for  deposition,  , 

fdep  '  r(EtfEa>  x  p  <416> 


where  p  is  a  random  number  uniformly  distributed  on  the  interval  (0, 1) 
and  r(E^pEa)  is  the  total  age  decrement  as  used  in  Eq.  (4.11).  Cor¬ 
responding  to  this  age  decrement  is  a  slowing  down  time  t^r  so 

that  we  may  determine  the  total  time  for  energy  deposition 


tj  +  t 


dep 


'4.  17) 


Furthermore,  the  spatial  position  for  energy  deposition  is 
determined  by  sampling  a  distance  r  from  Eq.  (4.  11)  with  r  =  r 

and  using  Eq.  (4. 12)  to  find  the  final  spatial  coordinates. 


The  amount  of  energy  deposited  varies  depending  on  whether 
the  neutron  is  in  the  air  or  in  the  ground:  In  the  air,  energy  is  deposited 
from  neutron  absorption  and  elastic  scattering  processes.  The  neutron 
absorption  is  due  primarily  to  the  (n,  p)  reaction  in  nitrogen  with  the 
proton  energy  deposited  locally  at  the  absorption  site.  The  (n,  y  )  reaction 
is  also  considered  but  the  gamma  ray  energy  is  transported  away  from 
the  absorption  site  as  will  be  discussed  later.  The  energy  deposited 
in  the  air  due  to  absorption  is  scored  in  the  analytical  routines  as 

Edep  =  (Eb  +  Q)  X  FRNPX  Wb  d-P(  W )  (*• 18> 


where  E^  is  the  neutron  energy  before  collision,  Q  is  the  energy  released 
during  the  (n,  p)  reaction  in  nitrogen  (.  6264  MeV),  FRNP  is  the  fraction 
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of  the  absorption  cross  section  which  results  in  proton  emission  (~.  96). 
wb  is  the  neutron  weight  before  collision,  and  p(E^-Ea)  is  the  resonance 
escape  probability  for  the  neutron  in  aging  from  E^  to  Ea- 

There  is  also  energy  deposited  in  the  air  due  to  elastic  scatter¬ 
ing  reactions.  In  this  case,  we  have  determined  empirically  that  approxi¬ 
mately  half  of  the  neutron  energy  is  transmitted  to  electrons  in  this 
energy  range.  Therefore,  elastic  energy  deposition  in  the  air  is  scored 
in  the  analytical  routines  as 

Edep  *  <Eb  '  Ea> x  Wb  x  P<  V  V  x  FE  (4' 19) 

where  and  Ea  are  the  neutron  energies  before  and  after  aging, 
is  the  neutron  weight  before  collision,  p(E^-*Ea)  is  the  resonance 
escape  probability,  and  FE  represents  the  fraction  of  energy  converted 
to  electrons  or  .  50  in  this  case. 

In  the  ground,  the  only  neutron  energy  deposited  is  due  to 
elastic  scattering  reactions  with  absorption  being  neglected  in  this 
energy  range.  For  ground,  we  assume  full  conversion  of  the  neutron 
energy  to  electrons,  so  that  Eq.  (4.  IP)  is  valid  with  FE  =  1.  0. 

4.  3.  3  Thermal  Neutron  Transport 

When  a  neutron  has  reached  energies  below  the  thermal  group 
cutoff  energy,  E^.,  a  new  time-dependent  one-group  diffusion  model  is 
utilized  whenever  infinite  medium  diffusion  theory  is  valid.  The  new 
model  is  applied  at  spatial  positions  which  are  distant  from  the  air- 
ground  interface.  At  altitudes  or  depths  which  are  closer  to  the  inter¬ 
face,  a  one- velocity  Monte  Carlo  calculation  is  invoked  so  that  the  two- 
region  geometrical  complications  can  be  accurately  treated.  The  new 
coupled  Monte  Carlo- diffusion  theory  scheme  does  enable  tracking  of 
thermal  neutrons  back  and  forth  from  ground  to  air.  In  this  manner, 
the  new  model  is  a  significant  improvement  over  the  ground  thermal 


capture  treatment  of  last  year  in  which  neutrons  leaving  the  ground  were 
assumed  lost  from  the  system  and  were  no  longer  followed. 

The  validity  of  the  analytical  model  can  be  determined  on  the 
basis  of  solid  angle  and  thermal  diffusion  theory  arguments.  We  have 
found  that  thermal  neutrons  located  at  altitudes  in  air,  or  at  ground 
depths,  more  than  three  diffusion  lengths  from  the  air-ground  interface 
can  be  treated  by  infinite  medium  diffusion  theory.  Only  a  small  per¬ 
centage  of  those  neutrons  located  three  or  more  diffusion  lengths  from 
the  air-ground  interface  will  eventually  cross  the  interface  and  enter 
the  other  medium.  Therefore,  only  in  rare  instances  do  we  need  to 
reject  the  direction  chosen  for  a  particle  which  is  analyzed  by  diffusion 
theory.  The  bias  introduced  into  the  calculations  by  this  direction 
selection  process  should  be  negligible. 

Those  neutrons  which  are  less  than  three  diffusion  lengths  from 
the  interface  are  followed  with  one  velocity  Monte  Carlo  and  may  travel 
from  one  medium  to  the  other.  The  Monte  Carlo  treatment  is  continued 
until  the  particle  is  either  killed  on  the  basis  of  weight  criteria  or  until 
it  moves  to  a  position  more  than  three  diffusion  lengths  from  the  inter¬ 
face  where  the  analytical  model  may  be  applied.  The  Monte  Carlo 
calculations  are  generally  much  slower  than  the  analytic  routines, 
especially  in  the  ground  where  several  hundred  collisions  may  be 
required  before  the  neutron  can  be  disposed  of  on  weight  criteria. 
Fortunately,  however,  the  neutrons  do  migrate  on  each  Monte  Carlo 
collision  so  that  they  may  diffuse  out  to  distances  where  the  more 
efficient  analytic  routines  may  be  employed. 

The  infinite  medium  thermal  neutron  flux  distribution,  away 
from  the  air-ground  interface,  is  obtained  by  solving  the  time  dependent 


174 


STirrw  — i 


riu'n  M— niiliMt 


diffusion  equation  for  a  spherically  symmetric  point  source  Qq  at 
time  t  =  0, 

I  *£  = 

v  3t 

where  </>  =  0(r,t)  is  the  thermal  neutron  flux  and  the  cross  sections 
and  speed  v  in  Eq.(4.  20)  are  thermal  values. 

The  solution  to  Eq.  (4. 20)  is  given  by 


I  <f>  +  Q 

r  art  a 


5(r) 

0  4^ 


8  (t) 


(4. 


0(M) 


-I  vt  -(r2/4Dvt) 
Qq  e  a  v  e 

“ 


(4. 


where  r  is  the  distance  from  the  source  point. 

The  coupled  thermal  diffusion  theory-Monte  Carlo  method  is 
employed  in  the  following  manner.  A  test  is  made  initially  to  deter¬ 
mine  whether  a  thermal  neutron  is  within  three  diffusion  lengths  of  the 
air-ground  interface.  If  the  particle  is  within  three  diffusion  lengths, 
the  one  velocity  Monte  Carlo  method  is  used  to  track  the  neutron  which 
may  be  either  in  the  air  or  in  the  ground.  The  Monte  Carlo  method  is 
used  until  the  particle  is  either  killed  by  weight  criteria  or  until  the 
neutron  is  transported  far  enough  from  the  interface  so  that  the  analyt¬ 
ical  treatment  can  be  employed.  Once  the  neutron  is  in  the  diffusion 
theory  region,  the  time  and  distance  traveled  before  its  absorption  and 
removal  from  the  system  are  computed  by  sampling  from  Eq.  (4.21). 
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Before  performing  the  sampling,  the  neutron  flux  is  first 
normalized  so  that  it  becomes  a  true  probability  density  function. 

This  is  accomplished  by  setting 

4_r ^ 

*M)  =  Q -  1  <P(  r,t)  (4.22) 

where  0(r,  t)  is  given  by  Eq.  (4.  21).  We  then  have 


He 


dr  f (r,  t)  =1.0 


(4.  23) 


Since  there  are  two  variables,  r  and  t,  which  are  to  be  sam¬ 
pled  from  Eq.  (4.  22),  we  first  form  a  marginal  density  f^t)  and  a  marginal 
distribution  F^t),  from  which  we  sample  a  time  t.  These  functions 
are  defined  as  , 


f  Vt 

f  i  (t) s  /  £M)dr  =  1  ve  a 
'o  a 


(4.  24) 


'j(t)  =/*  fj(t’)dt' 

Jc\ 


=  1  -  e 


■I„vt 


(4.  25) 


After  we  have  sampled  a  time  t  from  Eq.  (4.  25)  we  form  the 
conditional  distribution  F2(r|t),  where 

F2(rlt)-/'^fl  dr'  . 

Jo  Jo 


(4.  26) 
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We  then  sample  a  distance  r  from  the  distribution  of  Eq.  (4.  26) 

for  a  given  time  t.  The  sampling  of  the  distance  r  is  accomplished 

(  7 ) 

by  the  rejection  method'  which  is  very  fast  and  has  an  efficiency  of 
about  40%  for  this  particular  distribution  function. 

The  distance  r  then  represents  the  total  distance  traveled 
by  the  neutron  from  its  initial  coordinates  in  the  diffusion  theory  region. 

,  Zg),  then  the  final  location  of 
the  neutron,  or  its  absorption  coordinates,  are  given  by  (x.,y.,z.), 
where 

x4  =  Xg  +  a'  •  r 

y4  =  y3  +  0'  *  r  (4. 27) 

z4  =  Zg  +  y'  ■  r 

and  a’,  p',  and  y’  are  the  components  of  an  isotropically  generated  unit 
vector,  and  r  is  the  distance  sampled  above.  In  the  unusual  event  that 
the  particle  crosses  the  air-ground  interface,  a  new  isotropic  vector  is 
chosen  and  the  coordinates  (x^>y^>z^)  are  recalculated. 

In  a  similar  fashion,  the  initial  time  of  appearance  of  the  neu¬ 
tron,  tg,  in  the  analytical  region  is  added  to  the  time  t  which  is  sam¬ 
pled  from  Eq.  (4.  25)  to  obtain  the  total  diffusion  time,  before  the 
particle  is  absorbed, 

t4  =  tg  +  t  (4.28) 


If  these  initial  coordinates  are  (Xg,yg 
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When  the  neutron  is  being  transported  in  the  one-velocity  Monte 

Carlo  region,  its  weight  is  reduced  on  each  collision  by  the  factor 

p  =  2  /l.  The  fraction  of  weight  lost  on  each  collision  is  incorporated 
s 

into  gamma  ray  weight.  Upon  absorption  of  the  neutron  in  the  diffusion 
region  part  of  the  neutron  weight  remaining  is  carried  off  by  gamma  rays. 

As  in  the  case  of  the  epithermal  energy  region,  neutron  energy  is 
deposited  in  the  thermal  group  in  air  due  to  the  (n,  p)  reaction  in 
nitrogen.  Elastic  scattering  energy  deposition  is  neglected  for  the  thermal 
group  since  it  is  extremely  small.  The  thermal  neutron  energy  deposited 
in  air  is  given  by 


Edep  =  Q  X  FRNP  X  Wb  X  (1_p)  (4.29) 

where  Q  is  energy  released  due  to  the  (n,  p)  reaction,  FRNP  is  the 
fraction  of  the  absorption  cross  section  in  air  which  results  in  proton 
emission,  wb  is  the  neutron  weight  before  collision,  and  p  is  the  non¬ 
absorption  probability.  Equation  (4.  29)  is  valid  for  the  Monte  Carlo 
region  and  is  also  true  for  the  diffusion  region  if  it  is  understood  that 
in  the  analytic  region  p-3. 

For  the  ground,  the  neutron  energy  deposition  is  neglected 
for  the  thermal  group  since  it  is  very  small.  Gamma  ray  energy 
deposition  for  both  ground  and  air  is  discussed  in  the  next  section. 
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4.  3.  4  Gamma  Ray  Transport 


In  this  section  we  discuss  briefly  the  new  analytical  techniques 
used  to  treat  gamma  ray  transport  following  thermal  neutron  absorption. 
There  are,  of  course,  also  gamma  rays  produced  by  fast  neutron  re¬ 
actions  such  as  inelastic  scattering,  charged  particle  processes,  etc. 
These  gamma  rays  are  dealt  with  using  the  techniques  described  in 
earlier  reports  ^  4  knd  will  not  be  discussed  further  here.  We  will 
concern  ourselves  only  with  gamma  rays  which  are  produced  by  neutron 
absorptions  in  the  thermal  group.  Gamma  rays  which  are  produced  by 
epithermal  neutron  absorptions  are  tracked  by  Monte  Carlo  routines  for 
convenience.  The  analytical  model  described  below  could  in  principle 
be  extended  to  treat  these  photons  in  a  straightforward  manner. 

The  method  used  to  treat  gamma  rays  produced  by  thermal 
neutron  absorptions  is  bai  ed  on  combining  an  analytical  buildup  factor 
model  with  the  conventional  Monte  Carlo  method  used  in  the  past.  The 
buildup  factor  model  can  be  applied  to  gamma  ravs  which  originate  at 
absorption  sites  which  are  distant  from  the  air -ground  interface.  For 
gamma  rays  which  are  born  at  positions  closer  to  the  interface,  the 
previous  Monte  Carlo  techniques  are  employed.  In  this  regard,  the 
gamma  ray  transport  from  thermal  neutron  reactions  is  treated  in  a 
fashion  similar  to  the  epithermal  and  thermal  neutron  transport,  i.  e.  , 
analytical  models  are  used  away  from  the  air-ground  interface,  and 
Monte  Carlo  methods  are  applied  near  the  interface. 

Buildup  factor  approximations  have  been  applied  successfully 
to  many  gamma  ray  transport  problems.  These  approximations  are 
based  on  the  observation  that  the  calculation  of  the  unscattered  gamma 
ray  flux  is  usually  a  simple  matter  and  involves  only  an  exponential 
kernel.  A  buildup  factor  is  then  introduced  as  a  multiplicative  factor 


to  incorporate  the  effects  of  scattered  photons.  In  a  uniform  medium 
the  gamma  ray  flux  at  radial  position  r  due  to  a  unit  point  source  at 
the  origin  of  gamma  rays  of  energy  E  is  then  given  by 

-£(E)r 

*(r,E)  =  B0(r,E)  (4.30) 

Anr 

where  X(E)  is  the  macroscopic  absorption  coefficient  for  gamma  rays 
of  energy  E.  B«^(r,  E)  is  the  flux  buildup  factor  for  the  gamma  rays  of 
energy  E,  and  also  depends  on  the  distance  r  from  the  source  of 
photons . 

In  a  similar  manner  we  may  write  the  photon  current  at  r  due 
to  gamma  rays  of  energy  E  as 

-£(E)r 

J(r,E)  =  Bj(r,E)  ? - (4.31) 

4rr 

where  Bj(r,  E)  is  now  the  photon  current  buildup  factor  for  gamma  rays 
of  energy  E. 

We  are  ultimately  interested  in  determining  the  gamma 
ray  ionization  rate  as  well  as  the  various  Compton  currents.  The  ion¬ 
ization  rate  is  obtained  from  Eq.  (4.  30)  by  introducing  the  fluence  to  ion¬ 
ization  conversion  factor,  Ky(E), 

Ky(E)  =  E  <ra  +  E  apE  +  (E-1.02)  app.  (4.32) 

where  K  y  (E)  is  discussed  in  detail  in  Refs.  1  and  3. 


180 


Using  Eq.  (4.  32)  in  Eq.  (4.  30)  we  obtain  the  gamma  ionization  rate 
D(r)  =  J  Ky  (E)  0  (r,  E)  dE  (4.33) 

where  it  is,  of  course,  also  implied  that  <£(r,E)  and  D(r)  are  time 
dependent  quantities . 

The  radial,  or  horizontal,  Compton  currents  can  be  obtained 
from  Eq.  (4.  31)  by  introducing  the  mean  forward  electron  range  R,(E) 

I  0  1 

and  integrating  over  energy'  d/ 

-Rf(E) 

Jc(r)  =  e7  rnr  J(r>E)dE  <4-34) 

c'  ' 

where  A  (E)  is  the  Compton  mean  free  path  length. 

v 

The  currents  as  determined  by  Eq.  (4.  34)  are  defined  with  respect 
to  the  birth  site  of  the  photon.  To  determine  the  true  radial  or  horizon¬ 
tal  components,  we  must  translate  these  results  back  to  the  origin  of 
the  coordinate  system.  When  this  is  done  the  true  radial  or  horizontal, 
as  well  as  the  theta  or  vertical,  components  are  easily  determined. 

Before  Eqs.  (4.  30)  and  (4.  31)  can  be  used  for  the  photon  flux  and  cur¬ 
rent,  expressions  must  be  derived  for  the  flux  and  current  buildup  fac¬ 
tors,  B<£(r,  E)  and  Bj(r,  E).  We  have  derived  polynomial  expressions 
for  B0(r,  E)  and  Bj(r,E)  which  should  give  reasonably  good  results. 

In  both  cases  we  have  assumed  that  the  buildup  factors  can  be  developed 
in  a  series  of  the  form, 

B(r,E)  =  1  +  a(E)  •  (2(E)-  r)  +  b(E)  •  (X(E>  r)2  (4.3b) 
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where  a(E)  and  b(E)  are  coefficients  to  be  determined  for  each  photon 
energy  of  interest,  and  (2(E).  r)  is  the  photon  distance  from  the  sou.rce 
in  units  of  the  mean  free  path. 

We  have  determined  the  coefficients  a(E)  and  b(E)  for  both  the 
flux  and  current  buildup  factors  for  a  number  of  photon  energies.  These 
coefficients  were  obtained  by  performing  least  squares  fits  to  Monte 
Carlo  gamma  transport  data  in  air  supplied  by  V.  W.  Pine  of  AFWL. 

The  AFWL  data  span  the  energy  range  from  0.  5  MeV  to  8.  0  MeV.  We 
have  extrapolated  coefficients  obtained  from  these  results  to  higher 
energies  since  photons  in  air  at  energies  as  high  as  11. 0  MeV  are  of 
interest. 


We  have  compared  our  polynomial  expressions  for  the  buildup 
factors  with  some  other  derived  resultsv  '  to  ensure  accuracy.  The 
buildup  factor  expressions  derived  from  the  air  transport  calculations 
have  also  been  used  for  the  ground  transport  calculations.  The  buildup 
factor  calculations  are  influenced  primarily  by  the  atomic  number  of 
the  transport  medium  so  that  differences  between  air  and  ground  should 
be  small. 

After  a  neutron  absorption  site  has  been  determined  from  the 
epithermal  and  thermal  neutron  transport  routines,  a  photon  energy  is 
selected  from  the  appropriate  capture  gamma  spectrum.  Capture  gam¬ 
ma  spectra  for  both  air  and  ground  are  taken  from  the  latest  ENDF/B 
tabulations.  Epithermal  capture  gamma  spectra  are  assumed  to  have 
the  same  shape  as  the  therm  d  spectra. 

Tests  are  then  made  to  decide  whether  the  analytical  or  Monte 
Carlo  gamma  ray  transport  routines  should  be  used.  In  our  calculations 
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we  have  used  the  analytical  buildup  factor  models  for  gamma  rays  which 
are  born  more  than  three  mean  free  paths  from  the  air-ground  interface. 
Photons  closer  to  the  interface  are  treated  using  the  previous  Monte 
Carlo  CHARGE  routine.  All  gamma  rays  which  are  formed  from  epi¬ 
thermal  neutron  absorption  are  also  followed  using  the  CHARGE  routine. 

For  photons  which  are  distant  from  the  interface  we  can  sample 
the  distance  r  traveled  by  the  photon  before  its  demise  from  normalized 
probability  density  functions  of  the  form 


<4(r,E)  *  4irr2  B(r,E)  — 


(4.  36) 


where  N  is  a  normalization  factor,  and  B(r,E)  is  either  a  flux  or  cur¬ 
rent  buildup  factor  as  given  by  Eq.  (4.  35).  In  actual  calculations,  we  have 

_yr 

chosen  to  sample  r  from  the  exponential  distribution  e  and  then 
modify  the  weight  by  N  x  B(r,  E)  for  both  flux  and  current.  In  this  way 
the  ionization  and  Compton  currents  are  both  scored  at  the  same  spatial 
position,  but  the  magnitudes  differ  because  of  the  different  weight 
modifications. 

After  the  distance  r  is  sampled,  we  obtain  the  final  coordinates 
(Xg,yg,Zg)  from  the  initial  coordinates  (x^,y^,z^)  by  generating  a  unit 
isotropic  vector  with  components  (a",/3’\  y”)  and  forming 

xc  =  x.  +  a"  •  r 

5  4 

y5  =  y4  +  0"  •  r  (4.  37) 

ze  =  z.  +  y"  .  r 

5  4 

Particles  which  cross  the  interface  are  again  not  allowed, 
but  these  events  are  very  rare. 
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The  current  coordinates  at  (Xg,  y^,  Zg)  are  relative  to  the 
photon  birth  site  at  (x^,  y^,  z^).  By  translating  these  coordinates  to 
the  origin  of  coordinates  for  the  system,  we  easily  can  obtain  the 
true  Compton  currents. 

It  is  also  possible  for  us  to  generate  K  photons  per  neutron 
absorption  and  modify  the  weight  of  each  by  dividing  by  K.  We  have 
found  that  K  =  10  gamma  rays  per  thermal  absorption  to  be  a  reason¬ 
able  number  which  improves  our  statistical  results  for  the  analytical 
photon  transport  calculations  without  undue  computational  effort. 

It  has  been  assumed  that  the  photon  travels  from  its  birth 
position  to  its  final  positlr,'  instantaneously.  This  is  a  good  approxi¬ 
mation  since  the  timt  *'oq<  .  ";d  for  neutron  absorption  is  orders  of 
magnitude  longer  than  the  photon  t  ransit  time. 

4.  4  BIASING  TLCHN1QWS 

There  are  several  Mi  nte  Carlo  biasing  techniques  utilized  in 
05RNIES  to  reduce  compui  r  running  times  and  to  obtain  results 
which  would  otherwise  br  difficult  to  attain  using  conventional  methods. 
The  amount  of  biasing  used  in  the  present  version  of  the  code  has  been 
reduced  somewhat  from  earlier  versions  and  has  been  replaced  in  part 
by  analytical  routines. 

However,  to  obtain  reasonable  statistical  results  at  great  dis¬ 
tances  from  the  source,  path  length  stretching  is  still  required.  For 
neutrons  in  air,  stretching  is  used  for  particles  of  energy  greater  than 
an  epithermal  cutoff  energy,  E  =  0.  Ill  MeV.  Path  le.gth  stretching 
is  also  used  for  gamma  rays  in  the  air  to  attain  distant  results.  A  more 
complete  discussion  of  these  techniques  is  given  in  Ref.  2. 

(3) 

The  angle  of  scattering  biasing  discussed  in  last  year's  work 
has  been  eliminated  since  we  are  now  interested  in  neutrons  moving 
back  and  forth  from  the  ground  to  the  air.  In  the  earlier  calculations 
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the  scattering  angle  biasing  eliminated  many  of  the  then  unwanted 
neutrons  which  were  traveling  away  from  the  ground.  Source  direc¬ 
tion  biasing,  as  discussed  previously,  has  been  retained,  but  the 
source  energy  biasing  has  been  eliminated  with  the  new  energy  band 
version  of  the  code. 

The  previous  limited  scoring  and  capture  gamma  ray  produc¬ 
tion  biasing  schemes  have  also  largely  been  eliminated  and  replaced 
with  new  techniques.  For  neutrons  in  the  epithermal  range, 

E  s  e  £  E, ,  neutron  absorption  is  simulated  by  multiplying  the  neu- 

Vs  l 

tron  weight  by  the  nonabsorption  probability  at  each  collision.  It  is 
not  efficient,  however,  to  produce  and  track  capture  gamma  rays  at 
each  of  these  numerous  collisions.  Therefore,  an  input  quantity, 
RREPI,  is  used  as  the  probability  of  producing  capture  gamma  rays 
in  both  the  air  and  in  the  ground.  These  photons  which  are  produced 
then  have  their  weights  adjusted  to  reflect  the  modified  production 
rate. 

For  neutrons  in  the  thermal  group,  E  ^  Et,  limited  capture 
gamma  production  is  again  utilized,  with  another  input  quantity, 
RRTHM,  reflecting  the  fraction  of  times  a  gamma  ray  is  produced 
per  thermal  neutron  collision.  The  energy  deposited  by  thermal 
neutrons  is  also  only  scored  a  fraction  RRTHM  of  the  time  to  avoid 
wasting  large  amounts  of  computer  time  on  tallying  small  quantities. 
The  weights  of  both  the  capture  gamma  rays  and  the  neutrons  are,  of 
course,  modified  to  reflect  the  limited  photon  production  and  neutron 
scattering. 

It  should  be  noted  that  the  limited  photon  production  and 
scoring  is  used  for  particles  which  are  being  tracked  with  the  time 
consuming  Monte  Carlo  techniques.  The  analytical  methods  do  not 
include  biasing  routines  used  at  present  since  they  are  substantially 
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more  efficient  than  the  Monte  Carlo  approach.  We  also  mention  that 
it  is  possible  to  vary  the  number  of  photons  produced  per  collision 
in  the  epithermal  and  thermal  ranges,  through  RREPI  and  RRTHM, 
as  well  as  in  the  fast  range  through  NPHOT.  In  this  way  we  may 
adjust  the  average  weight  carried  off  by  photons  in  each  of  the  three 
ranges  so  as  to  obtain  lower  statistic?!  variance. 
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NUCLEAR  CROSS  SECTIONS 


4.  5 

4.  5. 1  Introduction 

In  these  calculations  we  have  assumed  that  the  air  was  composed 
of  nitrogen  and  oxygen  with  atom  percentages  of  79  and  21  percent,  re¬ 
spectively.  The  air  density  was  taken  to  be  1.29  mg/cc.  The  ground  was 
assumed  to  be  composed  of  hydrogen  (16%),  oxygen  (57%),  aluminum  (8%), 
and  silicon  (19%),  with  a  density  of  1. 7  g/cc. 

The  data  utilized  in  our  present  work  are  based  entirely  on  the 
most  recent  ENDF/B  tabulations  for  both  neutron  and  gamma  ray  produc¬ 
tion  cross  sections.  Gamma  ray  interaction  cross  sections  are  those 
used  in  our  previous  studies  and  are  sufficiently  well  known  for  EMP 
applications.  Since  the  neutron  and  gamma  ray  production  data  for  pre¬ 
vious  EMP  source  calculations  were  based  on  earlier  data  sets,  we  dis¬ 
cuss  below  some  of  the  differences  in  the  new  and  old  cross  section 
tabulations.  Furthermore,  some  estimates  of  the  sensitivity  of  EMP 
source  calculations  to  changes  in  nuclear  cross  sections  are  also  pre¬ 
sented. 

Preparation  of  new  cross  section  sets  for  use  with  05RNIES  is 
a  fairly  time  consuming  task  requiring  several  man-weeks  of  effort.  A 

number  of  processing  codes  must  be  employed  before  ENDF/B  Library 
data  can  be  placed  in  formats  useable  by  05RNIES.  Initially  the  RIGEL 
code  was  used  to  convert  the  ENDF/B  data  tapes  from  BCD  modes  to 
binary  modes.  RIGEL  was  also  used  for  combining  data  from  several 
data  tapes  onto  a  single  tape.  The  SAD  code  was  used  to  convert  data 
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from  tabulated  angular  scattering  probability  distributions  to  Legendre 
polynomials  for  hydrogen.  Code  678  was  used  to  create  neutron  scatter¬ 
ing  cross  section  tapes  and  the  PHI  angular  distribution  tapes,  as  well 
as  preparing  the  master  cross  section  data  tapes.  The  QDE  and  AVKER 
codes  were  used  to  create  energy  deposition  data  for  air  and  ground  col¬ 
lisions,  respectively.  Several  other  codes  were  also  used  for  creating 
gamma  production  probability  distributions  for  discrete  and  continuum 
gamma  cross  sections. 

4.5.2  Specific  Nuclear  Cross  Sections 

A.  5. 2.1  Hydrogen 

The  new  data  for  hydrogen  comes  from  Modification  1  of  the 

ENDF/B-in  Library  (February  1972).  The  new  hydrogen  cross  sections 

(3) 

should  be  very  similar  to  those  employed  in  previous  EMP  studies v  '. 

The  new  data  differ  primarily  in  format  changes  made  to  enable  the 
new  results  to  conform  with  ENDF/B-IH  library  specifications,  and 
in  modifications  to  add  photon  interaction  data.  The  photon  production 
data  for  hydrogen  are  simple,  consisting  of  only  a  single  neutron  capture 
gamma  ray  emitted  at  approximately  2.23  MeV. 

A.  5. 2. 2  Nitrogen 

The  nitrogen  neutron  and  gamma  production  cross  sections  were 
taken  from  Modification  3  of  the  ENDF/B-III  tabulation.  The  latest 
modification  again  was  made  to  enable  the  current  DNA  library  to  con¬ 
form  with  ENDF/B-IH  Library  specifications.  The  new  cross  sections 
for  nitrogen  differ  primarily  in  the  inelastic  scattering  data.  The  total 
cross  sections  for  the  new  tabulation  are  in  reasonably  good  agreement 
with  the  older  results.  The  inelastic  cross  sections,  however,  are 
generally  a  factor  of  two  or  so  lower  than  the  older  data.  The  new  gam¬ 
ma  production  cross  sections,  which  include  the  effects  of  inelastic 
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neutron  cross  section  changes,  also  differ  markedly  from  the  earlier 
ENDF-UNC  data/  9  * 

The  new  tabulation  also  includes  detailed  listings  of  inelastic 
scattering  level  and  charged  particle  decay  schemes.  For  nitrogen  we 
are  considering  some  26  discrete  inelastic  scattering  levels  and  no 
inelastic  continuum.  The  first  12  inelastic  levels  decay  by  conventional 
means  with  emission  of  gamma  rays.  The  14  higher  levels,  above  9 
MeV,  decay  by  emission  of  protons  in  addition  to  gamma  rays.  For 
the  standard  inelastic  scattering  processes  the  energy  deposited  by  the 
recoiling  nitrogen  atom  is  given  by 

Edep  =  Em-Eout®)-Q'  <438) 


where  is  the  neutron  energy  before  collision  and  EQut  (Q)  is 

the  energy  after  collision  and  is  dependent  on  the  level  excitation  energy 

Q. 

For  the  cases  of  inelastic  scattering  followed  by  charged 
14  13 

particle  emission,  »N  (n,  n'pLC  ,  the  total  energy  deposited  locally, 

14  0 

by  the  recoiling  atom  and  the  proton,  is  given  by 


(Q)  -  BE, 


(4.  39) 


where  BE  here  is  the  binding  energy  of  relative  to  gC^-^  +  jP*. 

It  is  assumed  that  the  final  carbon  state  is  the  ground  state  since  no 
gamma  rays  from  excited  carbon- 13  states  are  in  the  inelastic  ENDF/B 
tabulation. 


In  addition  to  the  inelastic  neutron  scattering  levels  which 
decay  by  charged  particle  emission,  a  number  of  conventional  absorp¬ 
tion  and  charged  particle  emission  reactions  are  also  considered  in 
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14  14 

the  new  cross  section  set.  For  „N  (n,  p)ftC  processes,  cross 

14  0 

sections  are  given  for  one  final  gC  ground  state  process  and  for 
four  transitions  to  excited  carbon  levels  with  eventual  gamma  ray  decay. 
For  deuteron  emission  processes,  ^N^(n,  d)gC^,  one  final  gC^  ground 
state  and  three  excited  carbon- 13  states  are  considered.  The  excited 

states  again  decay  by  gamma  emission.  For  triton  decay  processes, 

14  12  12 

(n,  t)gC  ,  one  ground  state  and  one  excited  gC  state  are  treated. 

(n,  a  )  reactions  are  considered  with  a  final  ground  and  also  ten  excited 

gB**  states  which  then  decay  by  gamma  emission. 

Thermal  neutron  (n,  y  )  capture  processes  in  nitrogen  are  also 
being  considered  in  our  EMP  calculations  for  the  first  time.  These 
cross  sections  are  also  available  from  the  ENDF/B  file.  The  multi¬ 
plicity  is  found  to  be  2. 166  and  the  average  gamma  emission  energy  is 
about  5  MeV.  This  results  in  a  total  gamma  source  energy  of  10.  83 

MeV  which  agrees  with  the  binding  energy  available  due  to  neutron 
14 

absorption  in  .It  is  assumed  that  the  thermal  capture  gamma 
spectrum  can  be  applied  to  higher  energy  ranges  also  since  spectral 
data  at  higher  energies  are  sparse.  This  should  not  introduce  any 
appreciable  error  since  the  thermal  group  is  most  important  for  (n,  y) 
reactions  and  furthermore  the  total  gamma  energy  available  of  10.  83 
MeV  remains  approximately  constant  at  epithermal  neutron  energies. 

4.  5.  2.  3  Oxygen 

The  oxygen  data  were  taken  from  the  ENDF/B  tabulation  of 

August  1971.  The  new  oxygen  data  include  recent  gamma  production 

(  9 ) 

cross  sections  which  differ  from  the  older  UNC  data  used  in  our  pre¬ 
vious  EMP  studies.  The  new  neutron  inelastic  scattering  cross  sections 
are  also  generally  lower  than  the  earlier  ENDF/B  results  used  in 
EMP  calculations. 
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The  new  ENDF/B  tabulation  includes  considerably  more  detailed 
charged  particle  and  inelastic  scattering  decay  schemes.  In  our  present 
calculations  we  have  considered  some  28  discrete  neutron  inelastic  scat¬ 
tering  levels  with  level  excitations  up  to  15.  25  MeV.  No  inelastic  con¬ 
tinuum  is  required  for  the  new  tabulation.  Of  the  28  inelastic  levels, 
only  6  are  conventional  (n,  n'y)  levels,  while  the  other  22  levels  decay 
by  either  proton,  alpha  particle,  or  alpha  plus  gamma  emission.  As  in 
the  case  of  nitrogen,  it  becomes  important  to  correctly  treat  the  charged 
particle  emissions  from  neutron  inelastic  scattering  levels  since  the 
charged  particles  deposit  essentially  all  of  their  energies  locally,  whereas 
gamma  rays  may  deposit  their  energies  over  much  greater  distances. 

For  the  six  conventional  inelastic  scattering  levels,  the  energy 
deposition  of  the  recoiling  oxygen  atom  is  again  determined  from  4.38. 

Edep  =  Ein  -  Eout(Q)  -  Q-  (4-38) 


where  Ein  is  the  neutron  energy  before  the  inelastic  scattering  collision, 
Eout(Q)  is  the  energy  after  collision  and  is  dependent  on  the  level 
excitation  energy  Q. 

The  energy  deposition  for  the  eight  inelastic  levels  which  decay 
by  alpha  emission  to  the  ground  state  of  carbon- 12  is  given  by 


E ,  =  E.  -  E  , 

dep  in  out 


(Q)  -BELC12), 


(4.40) 


where  BE  is  the  binding  energy  of  oxygen- 16  relative  to  breakup  into 
carbon- 12  plus  an  alpha  particle. 
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There  are  also  eight  inelastic  levels  which  are  assumed  to  decay 

first  by  alpha  emission  to  the  first  excited  state  of  carbon-12  and  then 

12 

decay  by  emission  of  4. 439  MeV  gamma  rays  to  the  gC  ground  state. 

The  local  energy  deposition  for  these  levels  is  computed  from 

Edep  ■  Ein  *  Eout<Q>  -  BE'  <4-41> 

where  Ein  is  again  the  incident  neutron  energy,  EQut  is  the  neutron  ener¬ 
gy  after  the  inelastic  scattering  event  for  level  Q,  and  BE  is  the  binding 
energy  for  oxygen  breakup  into  an  alpha  particle  plus  carbon- 12  plus  a 
4.439  MeV  gamma  ray. 

In  addition,  four  inelastic  levels  are  assumed  to  decay  by  proton 
emission  to  the  ground  state  of  nitrogen-15.  No  neutrons  are  energetic 
enough  in  our  calculations  to  result  in  inelastic  scattering  events  which 
decay  by  proton  emission  to  the  excited  states  of  nitrogen-15  with  sub¬ 
sequent  gamma  ray  decay.  The  local  energy  deposition  for  proton  emis¬ 
sion  to  the  nitrogen-15  ground  state  is  determined  from 

Edep  *  Ein  -  Eout(Q>  -  BE(7n15)’  <4-42> 

where  Ein  and  Eout(Q)  are  the  neutron  energies  before  and  after  scatter¬ 
ing  for  level  Q,  and  BE(  is  the  binding  energy  for  oxygen- 16  breakup 

into  nitrogen- 15  plus  a  proton. 

In  addition  to  the  gamma  rays  produced  by  inelastic  neutron 
scrttering,  gamma  rays  along  with  charged  particles  are  also  emitted 
in  absorption  processes.  Both  (n,py)  and  (n,  a  y)  processes  are  con¬ 
sidered  in  our  calculations.  The  (n,  y)  capture  reaction  in  oxygen  is 
very  small  and  has  been  neglected  in  our  studies. 
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4.5. 2. 4  Aluminum 

The  new  aluminum  data  were  taken  from  the  ENDF/B-II  Modi¬ 
fication  1  tabulation  of  July  1971.  There  were  35  discrete  inelastic  levels 
considered  in  the  latest  calculations.  No  inelastic  continuum  was  used, 
but  the  (n,  2n)  reaction  was  considered.  Local  energy  deposition  through 
charged  particle  reactions  was  treated  by  using  tabulated  kerma  values 
as  was  done  last  year. 

A  gamma  ray  continuum  was  included  in  the  new  calculations  as 
well  as  gamma  ray  production  cross  sections  due  to  neutron  inelastic 
scattering  and  capture  reactions. 

4.  5.2.5  Silicon 

The  new  silicon  data  were  taken  from  a  recent  ENDF/B  tabula¬ 
tion  and  include  revisions  by  R.  Q.  Wright  (2-05-71)  and  R.  R.  Kinsey 
(1-15-72).  The  new  neutron  and  gamma  ray  production  data  are  based 
on  the  M.  K.  Drake  evaluation.  Since  our  previous  EMP  silicon  data 
were  also  based  on  this  evaluation,  we  anticipate  little  change  in  the  sil¬ 
icon  cross  section  set.  We  have  compared  our  latest  silicon  data  with 
earlier  data  to  attempt  to  confirm  this  suspicion.  We  have  found  that 
the  new  and  old  total  cross  sections  differ  generally  by  less  than  about 
10%  from  18. 0  MeV  down  to  thermal  energies.  The  elastic  cross  sec¬ 
tions  also  are  in  good  agreement  and  the  inelastic  data  are  also  within 
about  10%  down  to  about  6  MeV. 

For  silicon  we  consider  21  discrete  neutron  inelastic  scattering 
levels  plus  a  continuum.  In  addition,  (n,  2n)  and  (n,n’p)  reactions  and  a 
gamma  ray  continuum  are  treated  in  these  calculations.  Local  energy 
deposition  for  silicon,  as  for  all  ground  materials,  is  handled  using  the 
AVKER  code  which  calculates  neutron  kerma  factors  which  are  then  used 
to  compute  dose  rates. 
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Cross  Section  Sensitivity  Studies 


Since  there  were  significant  cross  section  changes  between  the 
new  ENDF/B  data  and  older  cross  section  sets  used  in  EMP  studies,  a 
simple  set  of  cross  section  sensitivity  calculations  was  performed. 

The  changes  in  nitrogen  and  oxygen  neutron  inelastic  scattering  and 
gamma  ray  production  cross  sections  were  generally  the  most  substantial, 
and  the  sensitivity  studies  were  therefore  concentrated  on  these  materials. 

Neutron  and  secondary  gamma  ray  transport  calculations  were 

performed  in  infinite  air  using  the  latest  ENDF/B  neutron  and  gamma 

ray  production  cross  sections  for  nitrogen  and  oxygen.  A  second  set 

of  calculations  was  also  performed  using  the  earlier  data  sets  utilized 

(3) 

in  previous  EMP  studies  v  7  in  order  to  assess  the  sensitivity  of  neutron 
and  gamma  ray  transport  to  cross  section  perturbations.  The  calcu¬ 
lations  were  carried  out  using  the  one  dimensional,  time -dependent 
discrete  ordinates  code  TDA.  The  TDA  code  was  used  because  it 
required  less  computer  time  for  a  one -dimensional  problem  than  the 
Monte  Carlo  method.  More  detailed  two-dimensional  Monte  Carlo 
calculations  could  then  be  performed  if  the  TDA  results  indicated  further 
more  detailed  studies  were  warranted. 

A  point  source  of  neutrons  with  a  typical  thermonuclear 
source  energy  distribution  was  used  and  scalar  fluxes,  ionization  rates 
and  radial  currents  were  computed  using  the  two  sets  of  cross  sections. 
Results  were  obtained  out  to  distances  of  3100  meters  in  air  and  to 
times  of  approximately  100  nsec.  Thermal  neutron  capture  reactions 
were  not  considered. 

Detailed  plots  at  480,  1000,  2000,  and  3100  meters  of  radial 
Compton  currents,  neutron  and  gamma  ray  ionization  rates,  and  total 
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ionization  rates  were  made  to  compare  the  new  ENDF/B  data  results 
with  the  older  cross  section  data.  In  Figs.  4. 2-4. 5,  these  time -dependent 
ionization  rates  and  Compton  currents  for  the  two  data  sets  are  shown 
at  a  distance  of  480  meters  from  the  source.  Similar  results  were 
obtained  at  the  other  spatial  positions. 

In  Fig.  4.2  the  neutron  ionization  rate  for  the  new  data 
appears  to  be  substantially  higher  than  the  Mder  data  results,  parti¬ 
cularly  at  early  times,  and  again  at  later  times.  This  reflects  in  part 
the  lower  neutron  inelastic  scattering  cross  sections  used  in  the  new 
data  sets  which  result  in  less  spectral  degradation  of  neutrons  and 
larger  neutron  ionization  rates.  The  gamma  ray  production  cross 
sections  in  the  new  data  sets  are  substantially  lower  than  the  earlier 
data,  and  these  changes  are  exemplified  in  factor  of  two  reductions  in 
gamma  ray  ionization  rates  at  early  times  as  illustrated  in  Fig.  4.3. 
Figure  4.4  shows  the  total  ionization  rate  as  a  function  of  real  time. 

The  total  ionization  rate  as  computed  using  the  o’der  data  is  substantially 
higher  at  early  times  when  only  gamma  rays  are  contributing.  However, 
with  the  arrival  of  the  neutron  wave  front,  the  ENDF/B  data  calcu¬ 
lations  become  higher  than  the  older  results.  The  radial  Compton 
current  results  shown  in  Fig.  4. 5  generally  follow  the  shapes  of  the 
gamma  ionization  rates  with  the  older  data  results  being  substantially 
higher  at  early  times. 

The  cross  section  sensitivity  studies  imply  that  factor  of  two 
differences  in  quantities  such  as  radial  Compton  currents  and  total 
ionization  rates  may  be  introduced  by  using  the  latest  ENDF/B  data  for 
oxygen  and  nitrogen.  Changes  in  EMP  field  calculations  due  to  cross 
section  changes  should  be  somewhat  smaller. 
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Neutron  Ionization  Rate  (MeV/m  -sec-n) 


Figure  4.  2.  Comparison  of  neutron  ionization  rate  calculations 
using  ENDF/B  and  older  data  at  480  meters. 
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Figure  4.  4.  Comparison  of  total  ionization  rate  calculations 
using  ENDF/B  and  older  data  at  480  meters. 


Time  (seconds) 


Figure  4.  5.  Comparison  of  radial  Compton  current  calculations 
using  ENDF/B  and  older  data  at  480  meters. 


Brief  consideration  of  the  effect  of  cross  section  changes 
on  the  EMP  drivers  in  the  ground  was  also  made.  The  ground  used 
in  these  studies  was  composed  of  hydrogen  (16%),  oxygen  (57%), 
silicon  (19%),  and  aluminum  (8%).  The  new  ENDF/B  data  for  hydrogen 
should  not  differ  from  the  data  used  in  earlier  EMP  calculations.  The 
oxygen  data  are  considerably  different  and  the  TDA  numerical  calcu¬ 
lations  in  air  discussed  above  imply  that  factor  of  two  maximum  changes 
in  ionization  rates  and  compton  currents  may  result  with  the  new  data. 
Differences  in  ground  induced  fields  due  to  changes  in  oxygen  data 
should  be  similar  to  those  found  in  the  infinite  air  calculations.  The 
silicon  and  aluminum  components  of  the  ground  have  new  ENDF/B 
neutron  and  gamma  ray  production  cross  section  sets.  In  order  to  deter¬ 
mine  the  effects  of  these  cross  section  changes  on  EMP  drivers,  com¬ 
plicated  and  time  consuming  two-dimensional  transport  calculations 
would  be  required.  Therefore,  we  decided  to  simply  compare  the 
old  and  new  data  for  silicon  to  see  if  further  transport  calculations 
would  be  needed.  The  aluminum  data  were  not  compared  in  detail  since 
A1  is  only  an  8%  constituent  and  should  not  affect  the  EMP  drivers  to 
any  great  extent. 

For  Si  both  the  older  EMP  gamma  production  data  and  the  new 
ENDF/B  data  are  based  on  M.K.  Drake's  evaluation.  Therefore, 
gamma  ray  production  data  should  not  be  altered  with  the  new  ENDF/B 
data.  The  new  ENDF/B  neutron  data  are  also  based  on  Drake's  evalu¬ 
ation,  whereas  the  EMP  data  were  drawn  from  various  sources.  We 
have  compared  the  total  cross  sections,  as  well  as  the  total  elastic  and 
nonelastic  cross  sections  for  the  two  sets  of  data.  It  appears  that  dif¬ 
ferences  between  new  ENDF/B  neutron  data  and  the  EMP  data  are  not 
large.  For  example,  the  total  cross  sections  for  the  two  sets  generally 
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deviate  from  one  another  by  'v  +  10%.  There  are  a  few  points  where  the 
disagreement  is  somewhat  larger.  The  total  elastic  cross  sections  also 
seem  to  be  within  -10%  of  one  another,  while  the  nonelastic  cross  sections 
down  to  about  6  MeV  are  also  in  good  agreement.  Below  6  MeV  there 
are  larger  discrepancies  for  the  nonelastic  cross  sections. 

In  summary,  it  appears  that  the  new  ENDF/B  and  older  EMP 
neutron  and  gamma  ray  production  data  are  in  good  agreement  for  sili¬ 
con.  The  -10%  discrepancies  between  these  data  sets  are  relatively 
insignificant  when  compared  to  the  deviations  in  the  oxygen  data  sets. 
Therefore,  it  appears  that  a  full  transport  calculation  to  accurately 
compute  changes  in  EMP  drivers  due  to  changes  in  ground  cross  section 
data  is  not  warranted  at  this  time. 
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4.  6  INSTRUCTIONS  FOR  USING  05RNIES 

4.  6. 1  Q5RNIES  Input  Instructions 

The  input  information  required  to  run  05RNIES  is  discussed  in 
this  section.  The  subroutines  which  actually  read  the  input  data  are  also 
included  to  increase  the  user’s  familiarity  with  the  05RNIES  code. 

Rather  than  explain  in  detail  all  of  the  possible  options  for  05R  sub¬ 
routines,  the  instructions  herein  describe  the  input  as  these  subroutines 
are  used  in  05RNIES.  The  input  subroutines  used  here  include  INPUT, 

INELIN,  GATR,  MNIT,  SOURCE,  and  SSTART.  Input  instructions  for 
the  auxiliary  folding  and  plotting  code  05NPT  are  presented  in  Section 
4.  7. 1.  In  addition  to  the  input  cards,  which  are  discussed  in  this  section, 
tape  units  must  be  mounted  for  each  run.  These  are  discussed  in 
Section  4.  6.  3. 
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The  input  cards  for  05RNIES  are  as  follows: 


SUBROUTINE  INPUT 

CARD  A:  TITLE 

FORMAT  (20A4) 

TITLE  (dimensioned  20):  Hollerith  identification  of  the 
printed  output. 

CARD  B:  NSTRT,  NMOST,  NITS,  NQUIT,  EBOT,  NTHRML,  THPEED 

FORMAT  (415,  E10.  5,  15,  E10.  5) 

NSTRT:  The  number  of  neutrons  with  which  to  start  off 
each  energy  band. 

NMOST:  =  NSTRT 

NITS:  =  1 

NQUIT:  =  The  number  of  batches  to  be  run. 

EBOT:  The  energy,  in  eV,  below  which  neutrons  are 
considered  to  be  in  the  thermal  group. 

NTHERML:  -  1,  One- velocity  thermal  group 

CARD  C:  MEDIA,  NCONT1,  NCONT2,  NSTAPE 

FORMAT  (415) 

MEDIA:  =  2,  The  total  number  of  media  appearing  in  the 
system 

NCONT1:  =  3,  The  logical  number  of  the  system  data  tape 
(neutron  collision  cross  section  data). 

NCONT2:  =  NCONT1 

NSTAPE:  =  0 

A  set  of  the  cards  whose  descriptions  follow  D,  E,  and  F  must 
be  furnished  for  each  medium.  They  are  to  be  loaded  in  the  order:  Dp  Ep 
Fp  Dg,  Eg,  Fg,  with  the  media  1  (air),  2  (ground)  in  the  same  order  as  they 
are  listed  on  the  system  data  tape. 


CARD  D: 


(LF1  (LM,  M),  LM  =  1,  L),  where  M  =  media,  and  L  - 
number  of  scatterers  in  media  M 

FORMAT  (815) 

LF1  -  0:  Pick  from  isotropic  scattering  distribution 
=  N:  Pick  from  the  distribution  function 

n 

f*(M)  =  Y  <PK  8(vk  -h) 

and  use  the  Nth  scatterer  appearing  on  the 
phi  tape. 

=  +N:  Pick  from  the  anisotropic  distribution 
POO  =  (1  +  3f.  *0/2  and  use  the  Nth  set  of 
fj’s  on  the  system  data  tape. 

CARD  E:  (ASSES  (LM,  M),  LM  =  1,  L),  where  M  and  L  are  as 

defined  above  for  Card  D 

FORMAT  (7E10.  5) 

This  card  (or  cards)  carries  the  atomic  masses,  in 
amu,  of  each  scatterer  of  a  given  medium.  An  inelastic 
scatterer  is  represented  by  a  negative  mass. 

CARD  F:  SLOTH(M),  SLOPS(M),  SLOFS(M) 

FORMAT  (3E10.  5) 

SLOTH(M):  The  thermal  mean  free  path  for  the  medium  M 
in  cm. 

SLOPS (M):  The  thermal  nonabsorption  probability  for  the 
medium  M. 

SLOFS(M):  0. 0 

CARD  G:  ESOUR,  UINP,  VINP,  WINP,  WTSTRT 

FORMAT  (5E10.4) 

Not  used  since  subroutine  SOURCE 
specifies  initial  neutron  energy  and 
direction. 

The  statistical  weight  to  be  assigned  to 
each  source  neutron. 
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ESOUR:  =  0.) 
UINP:  =  0.  ( 

VINP:  =  0.( 

WINP:  =  O.j 

WTSTRT:  =  1. , 


CARD  H: 


CARD  I : 


CARD  J: 


CARD  K: 


XSTRT,  YSTRT,  ZSTRT,  NMED,  NREG 
FORMAT  (3E10.4,  215) 

The  coordinates  of  the  initial  position 
of  all  neutrons 

NMED:  =  0 
NREG:  =  0 

DUMMY 
FORMAT  (020) 

DUMMY:  The  octal  representation  of  the  initial  random 
multiplier.  If  equal  to  zero,  the  random  num¬ 
ber  function  starts  with  an  initial  random 
multiplier. 

NSPLT,  NKILL,  MXREG,  MAXGP 
FORMAT  (415) 

NSPLT:  0,  no  neutron  splitting 

NKILL:  =1,  neutrons  subjected  to  Russian  roulette 
(when  thermal) 

MAXREG:  =1 

MAXGP:  The  number  of  energy  groups  for  the  appli¬ 
cation  of  weight  standards. 

(WTS  (I),  1=2,  MAXGP) 

FORMAT  (7E10.  5) 

WTS:  Energy  values,  in  eV,  used  to  divide  the  energy 
range  of  the  problem  into  groups  for  the  applica¬ 
tion  of  weight  standards.  There  will  be  a  total 
of  (MAXGP-1)  values  of  WTS,  listed  in  descending 
order.  The  first  group  weight  standards  then  will 
apply  to  neutrons  having  energies  greater  than 
WTS(l).  The  second  group  of  standards  will  apply 


XSTRT:  =  0.  ) 
YSTRT:  =  0.  } 
ZSTRT:  =  0.  ) 


to  neutrons  with  energies  between  WTS(1)  and  WTS(2), 
and  so  on.  The  last  group  of  weight  standards  will 
apply  to  ail  neutrons  having  energies  below  WTS(MAXGP-l). 
Card  L  is  omitted  if  the  number  of  energy  groups  is 
one  or  zero  or  if  NSPLT  and  NKILL  are  both  equal 


to  zero. 
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CARD  L:  NGP1,  NDG,  NGP2,  NRG1,  NDRG,  NRG2,  WTHIH1, 

WTLOW1,  WTAVE1 

FORMAT  (15,  15,  15,  15,  15,  15,  E10.5,  E10.5,  E10.  5) 

These  parameters  designate  particular 
groups  and  regions  within  which  the 
weight  standards  given  by  g,  h,  and  i, 
below,  are  to  be  applied.  From  energy 
group  NGP1  to  energy  group  NGP2,  in¬ 
clusive,  in  steps  of  NDG,  and  from 
region  NRG1  to  region  NRG2,  inclusive, 
in  steps  of  NDRG  the  weight  standards 
to  be  applied  are  the  values  given  by 
g,  h,  and  i.  If  NGP1  is  set  to  zero,  all 
groups  will  be  given  the  weights  specified 
by  g,  h,  and  i,  and  if  NRG1  is  set  to  zero 
all  regions  will  be  given  the  weights  speci¬ 
fied  by  g,  h,  and  i.  As  many  cards  MM” 
as  are  required  may  be  used.  The  end 
of  the  loading  of  all  weight  standards  is 
signaled  by  a  separate  card  "M"  with 
NGP1  =  -1. 

g.  WTHIH:  Neutrons  having  weights  above  this  value 

in  the  groups  and  regions  indicated  will 
be  split. 

h.  WTLOW:  Neutrons  having  weights  below  this  value 

in  the  groups  and  regions  indicated  will 
be  subjected  to  Russian  roulette. 

i.  WTAVE:  If  a  neutron  is  not  killed  by  Russian 

roulette  it  will  be  assigned  a  new  weight 
equal  to  WTAVE. 

The  operation  of  these  parameters  may  be  clarified  by  a 
glance  at  the  following  FORTRAN  DO-loop: 

DO  1  I  =  NGPl,  NGP2,  NDG 

DO  1  J  =  NRG1,  NRG2,  NDRG 

WTHIH(I,  J)  =  WTHIH 
WTLOW  (I,  J)  =  WTLOW 
WTAVE  (I,  J)  =  WTAVE 
1  CONTINUE 

If  neither  splitting  nor  Russian  roulette  are  permitted,  card  "M"  is  omitted. 

206 


CARD  M: 


NSOUR,  MFISTP 
FORMAT  (215) 

NSOUR:  = 0 
MFISTP:  =  0 

CARD  N:  NPTAPE,  NTYPE 

FORMAT  (215) 

NPTAPE:  =  4,  The  logical  number  of  the  phi  tape 
NTYPE:  -  0 

SUBROUTINE  JOMIN 

CARD  As:  MED,  Z 

FORMAT  (15,  E10.5) 

MED:  Medium  with  lower  bound  of  Z. 

Z:  Z-axis  boundary  value,  in  cm.  The  Z  value  on 

the  MED  =  0  card  specifies  the  maximum  Z  value. 

CARD  B:  Blank  (no  region  geometry) 

CARD  C:  XL,  XU,  YL,  YU 

FORMAT  (4E10.5) 


XL: 

Lower  X  boundary, 

in  cm. 

XU: 

Upper  X  boundary , 

in  cm. 

YL: 

Lower  Y  boundary , 

in  cm. 

YU: 

Upper  Y  boundary, 

in  cm. 

SUBROUTINE  INELIN 

In  05RNIES, ,  each  scattering  reaction  that  an  element  may 
undergo  is  counted  as  a  separate  "scatterer.  M  Thus,  for  example,  an  ele¬ 
ment  with  5  discrete  levels,  elastic  scattering,  continuum  inelastic 
scattering,  and  an  (n,  2n)  reaction  contributes  8  separate  scatterers. 
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Information  about  all  scatterers  which  have  at  least  one  nonelastic  scat¬ 
tering  reaction  is  provided  through  INELIN. 

CARD  A:  COMM,  NINLS 

FORMAT  (11A6,  A4,  110) 

COMM  (dimensioned  12):  Description  of  INELIN  input 
(70  characters  are  allowed). 

NINLS:  Number  of  nonelastic  scatterers. 


CARDS  B  through  F,  below,  are  required  for  each  nonelastic  scatterer  (I). 


The  order  is  B 

gninls- 


V 


GV  B2’ 


G 


2’ 


B 


NINLS’ 


CARD  B:  LABEL,  NR (I),  NU(I),  NAPP(I) 

FORMAT  (5A6,  315) 

LABEL  (dimensioned  5):  Description  of  scatterer  (30 
characters  are  allowed). 

NR  (I):  =  1,  discrete  level  scatterer 

=  0,  otherwise 

NU(I):  =  0,  discrete  level  scatterer 

=  1,  otherwise 

NAPP(I):  Number  of  media  that  this  scatterer  is  in 

CARD  C:  (NMEDE  (J,  I),  NELE  (J,  I),  J  =  1,  NAPP  (I)) 

FORMAT  (415) 

NMEDE:  Medium  that  this  scatterer  is  in 
NELE:  Scatterer  number  in  that  medium 

If  NR  =  1,  then 

CARD  D:  Q(l,  I) 

FORMAT  (E10.4) 

Q:  Energy  level  of  this  scatterer  (in  CMCS). 

CARD  E:  None 

CARD  F:  None 
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CARD  G: 


None 


If  NU  =  1, 
CARD  D: 


If  NLAW 
CARD  E: 


CARD  F: 


CARD  G: 


then 


ELAW  (1,  I),  NLAW  a,  I) 

FORMAT  (E10.4,  110) 

This  card  describes  the  distribution  for  secondary  neutrons. 

ELAW:  Energy  threshold  for  this  scatterer  (in  eV). 

NLAW:  Indicates  distribution  for  secondary  neutron. 

Usually  NLAW  =  3  to  indicate  that  sampling  is 
to  be  made  from  Maxwellian  distribution. 


1 


FM  (1,  I),  NENG  (I),  NAE  (I) 

FORMAT  (E10.4,  2110) 

FM:  Multiplicity  for  secondary  emission.  (Example: 

FM  =  1.  for  inelastic  scattering;  FM  =  2.  for 
(n,  2n);  FM  =  3.  for  (n,  3n). ) 

NENG:  Number  of  energy  groups  for  this  scatterer. 

NAE:  Index  specifying  array  location  where  secondary 

energy  probability  distributions  are  stored  for 
this  scatterer. 

There  are  NENG  sets  of  F  and  G  cards  (one  F  card  and  as 

many  G  cards  as  necessary  for  each  set). 

ENG  (K,  NAE  (I)),  NETAB1  (K,  NAE (I)) 

FORMAT  (E10.4,  110) 

ENG:  Upper  energy  boundary  for  following  secondary 

energy  probability  distribution,  in  eV. 

NETAB1:  Number  of  secondary  energy  values  for  following 
probability  distribution. 


(ETAB1  (L,  K,  NAE  (I)),  PTAB1  (L,  K,  NAE  (I)),  L  =  1, 
NETAB1) 

ETAB1:  Secondary  energy  (eV) 

PTAB1:  Cumulative  probabilities. 
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CARD  E:  FM  (1,  I) 

FORMAT  (E10. 4) 

FM:  Multiplicity  for  secondary  emission. 

CARD  F:  TA  (I),  TB  (I),  TC  (I) 

FORMAT  (3E10.4) 

Nuclear  temperature^ (MeV)  for  Maxwellian  is  expressed 
as  T  =  TA  +  TB  *  \/E  +  TC  *  E,  where  E  is  the  incident 
neutron  energy  (eV). 


SUBROUTINE  GATR 

DO  55  1=  1,  2 

READ  (5,  1)  (A(J,  I),  J=  1,  11) 

READ  (5,  1)  (B(J,  I),  J=  1,  11) 

READ  (5,  1)  (AA(J,  I),  J=  1,  11) 

READ  (5,  1)  (BB(J,  I),  J=  1,  11) 

1  FORMAT  (6E12.  5) 

55  CONTINUE 

A  and  B  are  the  coefficients  for  the  buildup  factor  for  energy  deposition. 

AA  and  BB  are  the  coefficients  for  the  buildup  factor  for  current. 

1=1  is  air 
I  =  2  is  ground 

The  energy  group  boundaries  for  the  above  data  are  set  in  the  data  state¬ 
ments  in  GATR  (arrays  EA  and  EG).  For  further  details  see  Section  4.  3. 4. 
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u 

SUBROUTINE  XINIT 

CARD  A:  TIMMAX,  EPICUT,  RREPI,  RRTHM 

FORMAT  (4 El 0.5) 

TIMMAX:  The  maximum  time  limit  after  which  neutrons  are 

no  longer  followed,  in  sec 

EPICUT:  Upper  limit  of  epithermal  range,  in  eV. 

RREPI:  Probability  specifying  how  often  a  gamma 

particle  will  be  produced  while  tracking 
epithermal  neutrons. 

RRTHM:  Probability  specifying  how  often  a  gamma 

particle  will  be  produced  while  tracking 
thermal  neutrons. 

CARD  B:  NXTAP,  NHIS 

FORMAT  (216) 

NXTAP:  =  14,  Gamma  production,  Q,  and  AE  data  file 

NHIS:  Number  of  neutron  histories  per  energy  band 

CARD  C:  NTIMA,  NR  ADA,  NCOSA,  NZMA,  IWANT 

FORMAT  (516) 

NTIMA:  Number  of  time  bins 

NRADA:  Number  of  radial  bins 

NCOSA:  Number  of  angular  bins 

NZMA:  Number  of  horizontal  bins 

IWANT:  =  1 ,  Energy  is  not  deposited  for  first  neutron 

collision 

=  0,  Normal  condition 
>1,  Computes  energy  deposition  from  first 
neutron  collision  only. 

CARD  D:  (TIMC  (I),  1=1,  NTIMA  +  1) 

FORMAT  (6E12.5) 

As  many  card  D's  as  necessary 

TIMC:  Time  bin  boundaries  in  increasing  value,  in 

microsecs. 
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CARD  E: 


CARD  F: 


CARD  G: 


CARD  H: 


. ; 


(RADC  (I),  1=  1,  NRADA  +  1) 

FORMAT  (6E12.5) 

As  many  card  E's  as  necessary 

RADC:  Radial  bin  boundaries  in  increasing  value,  in 

meters. 

(ZADC  (I),  1=1,  NZMA+  1) 

FORMAT  (6E12.5) 

As  many  card  F's  as  necessary 

ZADC:  Horizontal  bin  boundaries  in  increasing  value, 

in  cm.  The  ground-air  interface  has  value  0.  0. 

(ANGC  (I),  1=1,  NCOSA+  1) 

FORMAT  (6E12.5) 

As  many  card  G's  as  necessary 

ANGC:  Angular  bin  boundaries  (cos  0)  in  increasing 

value,  with  0  being  horizontal  to  the  ground 
and  1.  0  in  a  vertical  direction. 

NPHOT,  WTLIM,  A,  NTHERM 
FORMAT  (112,  2E12.5,  112) 

NPHOT:  Number  of  photons  produced  at  each  photon  pro¬ 
ducing  collision  by  high  energy  neutrons. 

WTLIM:  Weight  below  which  Russian  roulette  is  played 
in  ROUTIN  to  determine  if  photon  is  followed. 

A:  Height  (cm)  of  the  neutron  source  above  the 

ground. 

NTHERM:  Number  of  capture  y’s  produced  for  each  thermal 
neutron  at  its  history  termination. 
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•JSKBaww 


0 

SUBROUTINE 
CARD  A: 

CARD  B: 


CARD  C: 

CARD  D: 


CARD  E: 


SOURCE 

NEB 

FORMAT  (15) 

NEB:  Number  of  energy  groups 

(E  (I),  I  =  1,  NEB  +1) 

FORMAT  (7E10.  4) 

As  many  card  B’s  as  necessary 

E:  Energy  group  boundaries,  in  MeV,  in  decreasing 
magnitude 

NUMA 

FORMAT  (15) 

NUMA:  Number  of  source  angular  distribution  bins 

(RMU  (I),  I  =  1,  NUMA  +1) 

FORMAT  (7E10.  4) 

As  many  card  D’s  as  necessary 

RMU:  Angular  bin  boundaries  (cos  6 ),  in  increasing 

numerical  order.  RMU  (1)  =  -1.  0  (straight  down ) 
and  RMU  (NUMA  +  1)  =  1.  0  (straight  up) 

(ASTAR  (I),  I  =  1,  NUMA) 

FORMAT  (7E10.  4) 

As  many  card  E's  as  necessary. 

ASTAR:  The  biased  source  intensities  for  the  angular 

bins.  The  assumed  true  source  angular  distribution 
is  isotropic. 
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SUBROUTINE  SSTART 


CARD  A:  IRESLT,  NRPT,  NEBC,  NBAC,  IDR 

FORMAT  (515) 

IRESLT:  05RNIES  results  tape. . 

NRPT:  Number  of  batches  per  results  tape  reel. 

NEBC:  Number  of  energy  groups  completed  for  batch 

presently  working  on. 

NBAC:  Number  of  batches  already  completed. 

IDR:  Run  identification  number. 

Typical  input  information  for  an  05RNIES  production  run  is 
shown  in  Table  7.  Some  of  the  information  comes  from  external  data 
files. 
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Table  7 


TYPICAL  INPUT  FOR  AN  05RNIES  PRODUCTION  RUN 


05HNIES  - 

£  MR  1972  • 

BAND 

VERSION  -  HOH  s 

0.00M  • 

IDaO  HOBaOOOM 

1000  1000 

1 

15 

2 

’.38 

1  2. 

2E5 

2 

i 

3 

0 

-1 

•?. 

0 

0 

0 

0 

0  0 

0 

0 

0 

0 

0 

0 

0  0 

0 

0 

0 

0 

0 

0 

0  0 

0 

0 

0 

0 

0 

0 

0  0 

0 

0 

0 

0 

0 

0 

0  0 

0 

0 

0 

0 

0 

0 

0  0 

0 

o 

0 

0 

0 

0 

0 

0  0 

14, 

006 

16. 

000 

-14. 

008 

-14,006 

-14,008 

-1 4,008 

-14,008 

-14, 

008 

-10. 

00  U 

-14. 

008 

-14,006 

-14.00- 

-  1  i ,  r.  6  6 

-1  4, OOP 

-14. 

008 

-1«. 

*08 

-14. 

008 

-14,006 

—  1  4,006 

•  1  4  ,  <’  (’6 

-1  4,008 

-14. 

008 

-14. 

006 

-14. 

008 

-14,008 

—14,006 

-14,008 

-14,008 

-14. 

008 

-  16, 

000 

•  16, 

000 

-16.000 

•16,000 

-16,000 

-16,000 

-16. 

000 

-  1  6  • 

00  0 

-  1  6, 

000 

-16.000 

-16, 0  r  O 

•16.00  f. 

-16,000 

-16, 

000 

-  1  6  . 

10  0 

-16. 

ooo 

-16,000 

-16,000 

•16,000 

-16,000 

-16. 

000 

-  1  6. 

000 

•  16. 

ooo 

-16.000 

-16,000 

-16,000 

-  16,000 

-  1  6, 

000 

1.H31F3 

,8528 

0,0 

-a 

•  2 

-3 

•  5 

0 

0 

0  0 

0 

0 

0 

0 

0 

0 

0  0 

0 

0 

0 

n 

0 

0 

0  0 

0 

0 

0 

0 

0 

0 

0  0 

0 

0 

0 

0 

0 

0 

0  0 

0 

0 

0 

0 

0 

0 

0  0 

0 

0 

0 

0 

0 

0 

0  0 

0 

0 

0 

0 

0 

0 

0  0 

0 

0 

0 

0 

0 

0 

0  0 

0 

0 

0 

0 

0 

0 

0  0 

0 

0 

0 

0 

0 

0 

0  0 

0 

u 

0 

0 

1. 

008 

16, 

000 

26 

.98 

28,086 

-16.000 

-16,000 

-16,000 

*16, 

000 

•  16, 

ooo 

-16, 

000 

•16,000 

-16,000 

- 1 6, 000 

-16,000 

- 16, 

000 

•16, 

ooo 

•16, 

000 

-16,000 

-16,000 

-16,000 

-16,000 

-16, 

000 

-16, 

000 

-16, 

ooo 

•16,000 

■  16,000 

-16,000 

-16,000 

■16,i 

000 

•16, 

ooo 

•16, 

ooo 

•16,000 

-26,98 

-26.98 

-26,98 

•  26 

,98 

•  26 

.98 

•  26 

.98 

-26,98 

-26,98 

-26.98 

-26,98 

•  26 

.98 

•  26 

,98 

•  26 

.98 

•26,98 

-26,98 

•26,98 

•26,98 

•  26 

.98 

-26 

.98 

-26 

.98 

-26,96 

-26,98 

-26,98 

-26,98 

•26 

.98 

•26 

.98 

•  26 

.98 

-26,98 

-26,98 

•26,98 

-26,98 

-26 

.98 

-26 

.98 

-26 

,96 

-26,98 

•26,98 

-28,086 

-28,086 

•  26,1 

086 

•28, 

086 

-28. 

086 

-28,086 

-28,086 

-28,086 

-28,086 

-28,i 

086 

-28, 

086 

-28. 

086 

•28,086 

-28,086 

•28.086 

-28.086 

•  ?8,l 

086 

-28, 

086 

-28, 

086 

•28,086 

-28,086 

-28,086 

-28,086 

-28,1 

086 

1. 

rt 

740 

.9912 

A  A  A 

A 

0,0 

A 

A  i 

A 

OtO  0.0  0,0  0,0  1,0 

0,0  0,0  0,0  0  0 

200032376454277*7 $41  RANDOM  NO. 


0  1 

1 

2 

2,38 

1  1 

1 

0 

0 

0 

1 , 0E*20 

1 .0E-20 

o 

rvj 

a 

o 

• 

nj 

2  1 

A 

2 

0 

0 

0 

1.0E+20 

1.0E-2 

1  .Ofc-l 

-1 

0  0 

4  0 
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Table  7  (Continued) 


2  •  250. 

l  •  o, 

0  5000. F.2 

■5 000, £2  5000. E2  -5000, £2  5000, £2 

I N 1 1  AST IC  INPUT  -  fcNOF/fl  CROSS  SECTION  DATA 
NITROGEN  (N.2NI)  0  1  t 

1  3 

1.13130E7  1 

2.  5  1 

1.1806E7  2 


6,25000*00  5,00000-01 
1.1750E7  1 

6,26000*08  2.37159-01 
1.300OE7  6 

6,25000*00  1.00691-01 
5.13500*05  5,6/301-01 
1.500067  15 

6,25000*044  1  ,52501-02 
5,1350O*(j 5  2,05392-01 
1.19500*06  5,69675-01 

2,nno0o*06  «,n6gp7-ni 

2.75000  +  06  9,94455  3  —  0 1 
1.70006/  23 

6.25000+  044  6.6244  /  3-03 
5,13500+  05  9. 2760/-02 
1.19500*06  2,59538-01 
2.00000+06  5,162/9-01 
2,75000*06  7,35020-01 
3.39900*06  8.95690-01 
9,25000*06  9,77762-01 
5,00000*06  9. 99886-01 

1ST  level  nitrogen 
1  a 
2.283E6 

2ND  LEVEL  NITROGFN 
1  5 

3.917E6 

3RD  LEVEL  NITROGEN 
1  6 

4I.8U8E6 

9 T H  LEVEL  NITROGEN 
1  7 

5.059E6 

5TH  LEVEL  NITROGEN 
1  8 

5.62JF6' 

6Th  LEVEL  NITROGEN 

1  9 

5, 796E6 

7TM  LEVEL  NITROGEN 
1  10 

6 , 1 448E6 

8TH  LEVEL  NITROGEN 
1  11 

6.3944E6 

9T H  LEVfL  NITROGEN 
1  12 

6.991E6 


1.87500*05  1,00000*00 

1.87500*05  7,11539-01 

1  .8  /500*05  3.360/3-01 
7,50000*05  9,8/680-01 

1,8/500*05  5, 1 000 0- 02 
7,50''00*„5  3.29  39  3-0! 
1,50000+06  7.10926-01 
2.250  0  0*06  9.-4  305  3-0  1 
3.00000*06  9,99928-OJ 

1  .87500*05  2.25M6-02 
/  ,  5 0 0 f.  0 ♦  6 5  |  ,52  3<>  3-9  1 

1.50000+06  3.63785-01 
2.25000*06  5.92026-01 
3.00000*06  7,9/392-01 
3. 75000*06  9.31389-01 
44,50000*06  9,90  386-01 
5.18200*06  1,00000*00 
1  0 


1  0 


1  0 


1  0 


1  0 


1  0 


1  0 


1  0 


1  0 


HDBsOOOM 

HOBsOOOM 


-  2/72 
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3,28000*05  1,00000*00 

3,26000*05  6,15861-01 
9,85000*05  1,00000+00 

3.2»  oro*„5  1 , Of  937-0  1 
9 ,  .  5  0  r  0 ♦ 0  5  8,02878-01 
1  ,  75'iOij*  06  6,1  230?-ni 

2,50Of0*Of.  0,  /7j  78-01 

3.1  89011*06  1  ,00000*00 

3 . 2  6  0  0  0  ♦  0  5  44,  /  7  1  6b-  0  2 
9, -50004  i*5  1.86922-01 
1  ,75000*06  44, 30  /82-01 
2,50000*06  6,65688-01 
3,10900*06  6,06832-01 
8,00000*06  9,56388-01 
8,75000*06  9,97386-01 


Table  7  (Continued) 


6.873F6 

9TH  LEVEL  OXYGEN 
1  33  2  8 

7.050E6 

5TH  LEVEL  OXYGEN 
1  39  2  9 

8.B21E6 

6TH  LEVEL  OXYGEN 
1  37  2  10 

9.999E6 

7TH  LEVEL  OXYGEN 
l  38  2  11 

9.707F6 

8Th  v-i  o>  ':  >n 

'  i<->  2 

io.27«r6 

9TH  LEVEL  OXYGEN 
1  35  2  53 

lO.Sdltb 

1  o f h  level  oxvr.fN 

1  3h  2  19 

10.983E6 

I  1  TH  LE  VH  OXYGEN 

1  <10  2  IS 

1  0  ,QF'.U  h 

12Tm  LEVtl  OXYGEN 
l  91  2  16 

5 i  ,116?  6 

I  ?  T  m  LEVEL  OXYGEN 
1  2  17 

11.337E6 

19TH  LEVEL  OXYGEN 
1  "3  2  18 

I I  ,«07F6 

1 5T H  LEVEL  OXYGEN 
1  99  2  19 

U.599E6 

1 6T H  LEVEL  OXYGEN 
1  95  2  20 

U.979E6 

1 7TH  LEVEL  OXYGEN 
1  59  2  21 

12.329E6 

1 8TH  LEVEL  OXYGEN  * 
1  96  2  22 

12.900E6 

19TH  LEVEL  OXYGEN 
1  97  ?  23 

12.683F6 

20TH  LEVEL  OXYGEN 
1  98  2  29 

12.896E6 

21ST  LEVEL  OXYGEN 
1  55  2  25 

13.037E6 

22Nf)  LEVEL  OXYGEN 
1  99  2  26 

1 3, 3? IF  6 

23fiD  LEVEL  OXYGEN 

1  50  2  27 


1  0  2 

1  0  2 

i  0  2 

1  0  2 

1  0  2 

1  0  ? 

1  0  2 

1  0  2 

1  0  2 

1  0  2 

1  0  2 

1  0  2 

1  0  2 

I  0  2 

1  0  2 

1  0  2 

1  0  2 

1  0  2 

1  0  2 

1  0  2 
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MitfaHiiMaifitii 


Table  ?  (Continued) 


l 

13.675F6 

29TH  LFVEL  OXYGEN  I  0  2 

1  56  2  28 

13.958E6 

2STH  LEVEL  OXYGEN  1  0  2 

1  51  2  29 

J9.292E6 

?6T H  LEVEL  OXYGEN  10  2 

!  57  2  30 

19.525E6 

27TH  LEVEL  OXYGEN  I  0  2 

1  52  2  31 

19.809F6 

28TH  LEVEL  OXYGEN  1  0  2 

1  53  2  32 

15.163E6 

ALUMINUM  ( N 1 2N)  0  1  1 

2  33 

1.35950E7  1 

2.  u  2 

1.3772E7  2 

5.00000+09  5.00000-01  l.bOfOO+OS  l.ooooo+oo 
1.9500E7  9 

5.O00O0+09  1.652OC-01  1.50000  +  05  5.02000-01  2.50000+05  8,36800-01 
3.50000+05  1,00000+00 
1.5500E/  19 

5.00000+09  3.20500-02  1.50000+05  1,03750-01  2.50000+05  1.85150-01 

3.50000  +  05  2.73150-01  50000  +  05  3,6380-1-01  5.500  01 +r 5  u.  53600-01 

6.50000+05  5.(11800-01  7.50000+05  6. 27150-01  0 , 50000  +  05  7.081  00-01 

9,50000+05  7. 83300-01  1. 05(100  +  06  8.52850-01  i.  15000  +  06  9,16550-01 

1,25000+06  9,73350-01  1,35000+06  1.00000+00 

1.6500E7  23 

5.00000+00  1.53500-02  1,50000+05  5,02000-02  2,50000+05  9,22500-02 

3. 50000  +  05  1.39650-01  <(,50000  +  05  1,90250-01  5,50000+05  2,02650-01 

6,50000+05  2,96700-01  7,50000+05  3.52000-01  8,50000+05  0,08050-01 

9,50000+05  0,60350-01  1,05000+06  5.20350-0J  1,15000+06  5,75350-01 

1.25000+06  6,29050-01  1,55000+06  6.81300-01  1,95000+06  7,31550-0! 

1,55000+06  7.79250-01  1,65000+06  8,23900-01  1,75000+06  8,69900-01 

1,85000+06  9.01600-01  1,95000+06  9,39250-0!  2,05000+06  9,63300-01 

2.15000+06  9,88950-01  2,25000+06  1.00000+00 

1ST  LEVEL  ALUMINUM  1  0  1 

?  39 

,  835E6 

2N0  LEVEL  ALUMINUM  1  0  1 

2  35 

1.007E6 

3RD  LEVEL  ALUMINUM  l  0  1 

2  36 

2.196E6 

9TH  LEVEL  ALUMINUM  1  0  1 

2  37 

2.709E6 

5TH  LEVEL  ALUMINUM  1  o  1 

?  38 

2.958E6 

6TH  LEVEL  ALUMINUM  1  o  1 

2  39 

2.977E6 
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Table  7  (Continued) 


7TH  LEVEL  ALUMINUM  1  0 

2  40 

3.60BE6 

6TH  LEVEL  ALUMINUM  1  0 

?  oi 
3.920E6 

9TH  LEVEL  ALUMINUM  .  1  0 

2  92 

4.029E6 

IOTh  level  ALUMINUM  l  0 

?  43 

9.379F6 

11TH  LEVEL  ALUMINUM  1  0 

2  04 

4. 4B3F6 

12'.H  LEVEL  ALUMINUM  1  0 

2  4b 
0  ,S56E6 

13TH  LEVEL  ALUMINUM  1  0 

2  46 

4.756E6 

IOTh  LEVEL  aluminum  i  o 

2  47 

S.J91E6 

1  ST H  LEVEL  aluminum  I  0 

2  OH 
S.69NE6 

1 61  m  LEVEL  AIUMIM.im  j  o 

2  04 

E.20HE6 

L7TH  LEVEL  ALUMINUM  t  0 

2  SO 
6.680F6 

18TH  LEVEL  ALUMINUM  J  (I 

2  51 

7.188E6 

19TH  LEVEL  ALUMINUM  1  0 

2  52 

7.696E6 

20TH  LEVEL  ALUMINUM  1  0 

2  53 

8.205E6 

2 1  ST  LEVEL  ALUMINUM  I  0 

2  54 

8.713E6 

22NI>  LEVEL  ALUMINUM  l  0 

2  55 

9.1ESE6 

23KD  LEVEL  ALUMINUM  1  0 

2  56 

9.656F6 

24TH  LEVEL  ALUMINUM  1  o 

2  57 

1 0 , J  6SE6 

2STH  LEVEL  ALUMINUM  |  o 

2  58 

10.673E6 

26TH  LEVEL  ALUMINUM  I  0 

2  59 
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Table  1  (Continued) 


1 1.109E6 


27Tw  LEVEL 
2  60 
11.617E6 

aluminum 

1 

0 

28TH  LEVEL 
2  61 
12.125F6 

aluminum 

i 

0 

29TH  level 
2  62 
12.634E6 

aluminum 

1 

0 

30TH  LEVEL 

2  63 

13.142E6 

ALUMINUM 

1 

0 

3 1 T H  LEVEL 

2  64 

13.650E6 

ALUMINUM 

1 

0 

32T H  LEVEL 

2  6b 

14.159E6 

ALUMINUM 

1 

0 

33T m  LEVEL 

2  66 
14.667E6 

ALUMINUM 

1 

0 

34 TH  LEVEL 

2  6/ 
lb ,  1  7T>F6 

ALUMINUM 

1 

0 

3bT H  IEVLL 

2  68 
tb.611E6 

Al  ^ l M i Ml J M 

1 

0 

SILICON  ( N « 
2  69 

6.7/0GE6 

2. 

-.602271  b 

?N) 

3 

,2lS60E3-.<42488E-2 

0 

1 

SILICON  (N 

2  70 

1.2000E7 

1 1 

-,?3H07E'6 

•  NPR I M  E  )  P 

3 

,29357E3-,13242E»1 

0 

1 

1ST  LEVEL 
2  71 

1.254E6 

SILICON 

1 

0 

260  LEVEL 
2  72 

1.763E6 

SILICON 

1 

0 

3RD  LEVEL 
2  73 

2.017E6 

SILICON 

1 

0 

4TH  LEVEL 
2  74 

2.217F6 

SILICON 

I 

0 

5TH  LEVEL 
2  75 

2.399E6 

SILICON 

1 

0 

6TH  LEVEL 
2  76 

3.035E6 

SILICON 

1 

0 

7TH  LEVEL 

SILICON 

1 

0 

2  77 

3.362E6 
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Table  7  (Continued) 

8TH  LEVEL 

SILICON 

l 

0  1 

2  78 

3.599E6 
9TH  LEVEL 

SILICON 

» 

0  1 

?  79 

3.726E6 

1  oth  level 

SILICON 

l 

0  1 

2  80 
3.744E6 
1ITH  LEVEL 

SILICON 

t 

0  1 

?  81 
4.580F6 
12TH  LEVEL 

SILICON 

l 

0  1 

2  82 
4.908E6 
13TH  LEVEL 

SILICON 

1 

0  1 

2  83 

6.217F6 
1«TH  LEVEL 

SIl  ICON 

i 

0  1 

?  84 

6.655F6 
15Th  level 

SILICON 

l 

0  1 

2  85 

6.835F6 

1  6 T H  LEVEL 

SILICON 

l 

0  1 

2  "6 
6.835E6 
17TM  LEVEL 

SILICON 

l 

0  1 

?  87 

7.5«UE  h 

1 8 T H  LEVEL 

SILICON 

l 

0  1 

?  88 

7,  J4«F6 

1 9 T H  LEVEL 

SILICON 

l 

0  1 

2  89 

7.743E6 
20TH  LEVEL 

SILICON 

l 

0  ! 

2  90 

7.889E6 

21  ST  LEVEL 

SILICON 

i 

0  1 

2  91 

8.216E6 

CONTINUUM  SILICON 

0 

1  1 

2  92 

4.0000E6 

1  • 

•  23272F  3 
.38415 

3 

.17802F3 

1,06546 

,6B193E.J 

,8645 

.7048 

,58017  .517 

4544 

.«1 

,36615 

.34 

.31346 

68507 

.15638 

,02877 

.00605 

1.362  -03-  .0012 

3,7702-03 

-  4,04 

-03-  4,3134-05- 

4.32  -0.3- 

4,330-03 

« 739 

,499 

,462 

,399 

.364  .369 

.355 

,340 

.326 

,326 

.326 

.0431 

,0365 

,00881 

,00638 

.004  ,0016 

-8,84  -04 

-2,07 

-03  -.00236 

-.00236 

-.00236 

.36415 

1,08346 

.8645 

.7048 

.58017  ,517 

4544 

.41 

,36615 

.34 

.31348 

68507 

.15638 

.02877 

,00605 

1,582  -03-  .0012 

3,7702-03 

•  4,04 

-03-  4,3134-03- 

4,32  -03- 

4.310-03 

AA  t 
AA  2 
BB  1 
BB2 
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Table  7  (Continued) 


1 734 

.499 

.462 

.  399 

,  .384 

,369  4 A  1 

.555 

,340 

.326 

.326 

,326 

AA  2 

.0431 

.0365 

,00881 

.00636 

1  .004 

,0016  HH  1 

•8.84  -oa 

•  2,07  -03 

-.00236 

-.00236  -.00236 

8B? 

1.0 

.11  If 6 

.1 

.1 

14  1000 

59  19 

5  9 

0 

HO9S000M 

.000 

.015 

.020 

.030 

.040 

,O6OHOfl*O0OM 

.080 

,100 

,140 

.200 

.300 

,  4  0  OHOHs  00  0M 

.600 

,600 

1.00 

2.00 

3.00 

4.O0HO6SO00M 

5.00 

7,50 

10,0 

20,0 

30.0 

5O.0HO4s000M 

75,0 

100. 

200. 

300, 

500, 

lOOO.HOrtsOOOM 

1500. 

3000. 

6000. 

l.M 

3.E4 

6.F  UHOBsOOOM 

1.E5 

3.F5 

6.E5 

1.F6 

HOBsOOOM 

00,0 

150. 

250. 

350, 

450. 

550 ,H0H»OOOM 

650. 

850. 

1050. 

1  300, 

1500, 

l7O0.HO«»000M 

1900, 

2100. 

2300. 

2700. 

3100. 

3500 .WOHsOOOM 

4000. 

5000, 

H(Jt3s000M 

•250, 

•  180. 

-120. 

-90. 

-50, 

-30 .HOHsOOOM 

-14. 

-6,0 

-2.0 

o.oo 

HODsOOOM 

.0 

.05 

.10 

.20 

.50 

1 ,0HO8s000M 

2  0 

1 ,£-10 

0,0 

10 

M03=OOOM 

13 

16.5 

15.0 

13.5 

g 

.05 

4.96 

2.35 

1.10 

.55 

.24 

.11 

.0316 

.0193 

00336 

.000023 

6 

KOHsOOOM 

•1.0 

0.00 

.05 

.10 

.20 

.50 

l.o  mfisooOM 

5.0 

1.0 

1.0 

1.0 

1.0 

1.0 

HOHsOOOM 

16  3 

0  12  0 

HOBsnoOM 

223 


■' — - — | f  ■ 


4.  6.  2  05RNIES  Output 


During  each  run  a  small  amount  of  printed  output  of  a  general 
nature  is  produced.  In  addition,  detailed  ionization  rate  and  Compton 
current  information  is  tabulated  and  plotted  on  microfilm.  The  inform¬ 
ation  printed  during  each  run  is  used  to  indicate  how  well  certain 
input  parameters  and  biasing  techniques  were  chosen.  Restart  inform¬ 
ation  is  also  printed  to  enable  the  user  to  restart  calculations  which 
were  aborted  for  any  reason. 

Figure  4.  6  shows  the  normal  output  produced  by  05RNDES. 

The  tracking  part  of  the  output  consists  of  the  total  weight  and  the 
weighted  average  position  for  all  neutrons  in  a  given  energy  band,  along 
with  the  last  random  multiplier,  after  all  neutrons  have  become  thermal, 
and  all  neutron  histories  have  been  terminated.  Also,  the  computer 
time  required  for  non-thermal  and  thermal  tracking  and  subsequent 
energy  deposition  is  noted.  At  the  end  of  each  energy  band  the  integer 
(ICOUNT)  and  floating  point  (PCOUNT)  counter  values  are  printed.  At 
the  end  of  a  batch  the  counter  values  summed  over  all  energy  bands  in 
that  batch  are  printed  as  KCOUNT  and  ZCOUNT  arrays  with  the  floating 
point  counters  being  normalized. 

The  KCOUNT  array  is  used  to  record  the  number  of  particles 
produced,  and  number  of  neutrons  whose  histories  are  terminated  in 
specific  ways,  and  the  number  of  collisions  of  a  given  type.  The  ZCOUNT 
array  records  total  energy  deposited  by  neutrons  and  gammas,  the 
weight  and  weighted  energies  of  neutrons  and  gamma  rays  and  various 
other  important  integral  quantities. 

Definitions  of  the  quantities  contained  in  the  KCOUNT  and 
ZCOUNT  arrays  are  listed  in  Tables  8  and  9. 
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Figure  4.6.  Typical  05RNIES  output. 


Table  8 


KCOUNT  Array 


Number 

Set  in 
Routine 

Quantity 

1 

BANKR 

Number  of  source  neutrons 

2 

BANKR 

Number  of  neutrons  escaping  from  system 

3 

BANKR 

Number  of  boundary  crossings 

4 

DOSEN 

Number  of  collisions  in  the  ground  analyzed 
by  DOSEN 

5 

BANKR 

Number  of  neutrons  whose  local  time  exceeds 
TMAX 

6 

05RNIES 

Number  of  neutrons  reaching  thermal  in  air 

7 

DOSEN 

Number  of  oxygen  inelastic  collisions  in 
the  air  analyzed 

8 

DOSEN 

Number  of  nitrogen  inelastic  collisions  in 
the  air  analyzed 

9 

DOSEN 

Number  of  continuum  (n,  n  )  nitrogen 
collisions  in  the  air  analyzed 

10 

DOSEN 

Number  of  (n,  2n)  nitrogen  collisions  in 
the  air  analyzed 

11 

DOSEN 

Number  of  elastic  nitrogen  collisions  in 
the  air  analyzed 

12 

DOSEN 

Number  of  elastic  oxygen  collisions  in  the 
air  analyzed 

13 

GATR 

Total  number  of  gammas  produced  by 
epithermal  and  thermal  neutron  reactions 

14 

NIES 

Number  of  collisions  analyzed 

15 

ROUT  IN 

Number  of  inelastic  gammas  produced 

16 

BANKR 

Number  of  neutrons  killed  in  05R  by 

Russian  Roulette 
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Table  ?  (Continued) 


Number 

Set  in 
Routine 

Quantity 

17 

BANKR 

Number  of  neutrons  surviving  Russian 
Roulette  in  05R 

18 

GATR 

Number  of  gammas  followed  by  GATR  for 
neutrons  terminated  by  RDIFF 

19 

05RNIES 

Number  of  neutrons  reaching  thermal  in 
ground 

20 

GATR 

Total  number  of  gammas  followed  by  GATR 

21 

BANKR 

Number  of  neutrons  taken  to  thermal  energy 
analytically 

22 

BANKR 

Number  of  neutrons  terminated  by  analytic 
techniques  in  RDIFF 

23 

BANKR 

Number  of  gammas  generated  by  thermal 
neutrons 

24 

BANKR 

Number  of  epithermal  collisions  (also 
number  of  non- analytic  epithermal  gammas) 

25 

BANKR 

Number  of  gammas  produced  by  analytic 
epithermal  neutron  reactions 

26  thru 

30 

Not  used 

KCOUNT  Checks 

KCOUNT  (1)  =  KCOUNT  (2+5+16+22) 

KCOUNT  (14)  =  KCOUNT  (4+7+8+9+10+11  +  12) 

KCOUNT  (1)  *  KCOUNT  (6  +  19) 

KCOUNT  (13)  =  KCOUNT  (24  +  25)  +  KCOUNT  (23)  +  N  THERM  *  KCOUNT  (22) 


NOTE:  For  conciseness,  KCOUNT  (6+19)  =  KCOUNT  (6)  +  KCOUNT  (19),  etc. 


Table  9 


ZCOUNT  Array 


Number 

Set  in 
Routine 

Quantity 

1 

BANKR 

DOSEN 

Energy  of  neutrons  outside  spatial  bins 

2 

DEP 

Energy  deposited  by  thermal  neutron  reactions 

3 

DOSEN 

Energy  deposited  by  fast  neutron  reactions, 

E  >E 

c 

4 

CHARGE 

Energy  deposited  by  gammas  produced 
by  fast  neutron  reactions,  E  >  Ec 

5 

ROUTIN 

Weight  of  inelastic  gammas 

6 

ROUT  IN 

Energy  of  inelastic  gammas 

-  7 

BANKR 

Weight  of  source  neutrons 

8 

BANKR 

Energy  of  source  neutrons 

9 

CHARGE 

GATR 

Energy  deposited  by  gammas  produced  by 
thermal  neutron  reactions,  E  <  E^. 

10 

GATR 

Energy  of  gammas  created  by  epithermal 
neutrons,  Ec?E>  E^. 

11 

GATR 

Energy  of  gammas  created  by  thermal  neutrons, 

ErE 

12 

DEP 

Part  of  2  deposited  in  the  ground 

13 

DOSEN 

Part  of  3  deposited  in  the  ground 

14 

CHARGE 

Part  of  4  deposited  in  the  ground 

15 

BANKR 

Energy  of  neutrons  whose  elapsed  local  times 
exceed  TMAX 

16 

CHARGE 

Energy  deposited  by  gammas  produced  by 
epithermal  neutron  reactions,  E£^  E^E^ 

17 

CHARGE 

Part  of  16  deposited  in  ground 

228 


Table  9  (Continued) 


Number 

Set  in 
Routine 

Quantity 

18 

05RNIES 

Weight  of  neutrons  in  air  which  reach 
thermal  energies 

19 

CHARGE 

GATR 

Part  of  9  deposited  in  ground 

20 

05RNIES 

Weight  of  neutrons  in  ground  which  reach 
thermal  energies 

21 

DOSEN 

Energy  deposited  by  epithermal  neutron 
reactions 

22 

23  thru 

30 

DOSEN 

Part  of  21  deposited  in  ground 

Not  used 

ZCOUNT  Checks 

ZCOUNT  (7)  =?  1 . 

ZCOUNT  (8)/ ZCOUNT  (7)  -  expected  energy  of  source  neutrons 
ZCOUNT  (4)  k  ZCOUNT  (6) 

ZCOUNT  (9)-  ZCOUNT  (11)  ^  .  04  *  10.  37  *  ZCOUNT  (18)  +  4.  923  * 
ZCOUNT  (20) 

ZCOUNT  (16)  -  ZCOUNT  (10) 

ZCOUNT  (2)  — .  624  *  ZCOUNT  (18)  *  .  96 
ZCOUNT  (12)  =  0. 

ZCOUNT  (3)  +  ZCOUNT  (6)  ~  ZCOUNT  (8) 
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4.  6.  3  Structure  of  the  Q5RNIES  Code 
4.  6.  3. 1  Flowcharts 

Figures  4.  7  -  4.  25  show  flowcharts  of  the  various  subroutines 
in  the  analysis  section  of  05RNIES.  Most  subroutines  in  the  05R 

(12) 

section  are  not  included,  as  they  have  been  documented  elsewhere. 
However,  a  few  05R  subroutines  either  are  not  included  in  Reference  12 
or  have  been  substantially  modified,  and  we  include  flowcharts  for 
these.  Among  these  are  BANKR,  SOURCE,  HINNY,  and  INELIN. 

Figure  4.  7  shows  a  block  diagram  of  what  we  refer  to  as  the 
05R  section  of  05RNIES.  This  section  of  the  program  generates  the 
neutron  histories.  Figure  4.  8  shows  a  block  diagram  of  the  analysis 
section  of  05RNIES;  it  performs  the  analysis  of  the  neutron  histories. 
Although  it  is  net  shown,  several  of  the  routines  in  the  05R  section 
are  also  used  by  the  analysis  section.  For  example,  CHARGE  which 
performs  a  complete  gamma  ray  transport  calculation  calls  the  sub¬ 
routine  GTISO,  which  is  also  called  by  several  routines  in  the  05R 
section.  The  remaining  pages  of  this  subsection  contain  flowcharts  of 
the  individual  analysis  subroutines  and  selected  05R  subroutines.  Each 
subroutine  flowchart  is  accompanied  by  descriptive  material  about  the 
subroutine. 

Except  for  BANKR,  the  05R  subroutines  are  almost  identical 
to  those  described  in  Reference  12.  We  have  made  some  minor  changes 
in  some  of  the  routines,  however.  These  routines  are: 

RWCON  -  Same  as  before,  except  it  now  rewinds 
NXTAP  and  sets  the  initial  values  of 
ETOPI  and  EBOTI. 
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05RNIES 


Figure  4.8.  Analysis  routines  in  05RNIES 


232 


READR  -  Same  as  before,  except  it  now  reads 
NXTAP  and  sets  ETOPI  and  EBOTI. 

NXTCOL  -  Neutrons  in  the  air  with  energies  above 

AIRCUT  have  their  path  lengths  stretched. 

A  number  of  subroutines  and  functions  are  called  only  by  sub¬ 
routine  KINNY.  We  do  not  include  flowcharts  for  each  of  these,  but 
below  is  a  brief  description  of  the  role  played  by  each  of  these  subroutines 
and  functions: 

AB  -  Obtains  ALPHA  and  BETA  for  a  given  target  nuclide. 
ALPHA  =  Mn/(Mn  +  Mt)  and  BETA  =  M,^(MN  +  MT),  where  MN 
is  the  mass  of  the  incident  neutron  and  is  the  mass  of  the  target 
nuclide. 

CHOS  -  Linear  interpolating  function.  Given  the  arrays  P 
and  Q  (of  size  NPQ),  where  Q  is  a  function  of  P,  and  P  is  in  decreasing 
numerical  order.  Then,  for  a  given  R  within  the  range  P(l)  to  P(NPQ), 
CHOS  will  determine  the  corresponding  Q  by  linear  interpolation. 

DCS  -  Computes  the  new  direction  cosines  for  the  neutron. 

ELAWS  -  Determines  the  speed  of  neutron  after  the  non¬ 
elastic  collision  exciting  continuum  levels. 

EVAP  —  Determines  the  neutron  energy  after  inelastic  scat¬ 
tering  by  evaporation  model. 

LEVELS  -  Determines  a  speed  used  in  calculating  the  new 
velocity  components  of  a  neutron  after  collision  with  an  inelastic  dis¬ 
crete  level  target  nuclide. 
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The  analysis  section  of  05RNIES  contains  several  additional 
subroutines,  for  which  we  do  not  include  flowcharts.  A  brief  descrip¬ 
tion  of  the  role  played  by  these  subroutines  is  as  follows: 

KLEIN  -  Routine  used,  to  select  scattering  angle  and  resulting 
energy  from  Compton  scattering  using  Klein-Nishina  cross  section. 

MOWER  -  Determines  the  smallest  power  of  2  greater  than 
or  equal  to  a  given  number. 

NINCS  -  Determines  which  value  in  the  NIES  cross  section 
data  array  is  to  be  used. 

XION  -  Computes  fraction  of  energy  deposition  in  N  and  O 
which  causes  ionization.  (From  fits  to  calculations  of  Lindhard, 
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Subroutine  Name: 


Function: 


BANKR  (NBNKID)  (Fig.  4.  9) 

Calls  subroutine  NIES  to  process  collision  information 
or  update  accumulative  results  arrays.  It  also  sets 
certain  ICOUNT  and  PCOUNT  values. 


Calling  Program:  05R,  MSOUR,  NTXCOL,  TESTW 

Subroutines  Called:  CHARGE,  DEP,  GATR,  NIES,  RAGE,  RDIFF, 

ROUTIN,  TAUM 

Input  Parameter:  NBNKID:  Specifies  what  action  BANKR  is  to  take. 

See  Ref.  12  for  details. 


Labeled  COMMON:  CALLE  ,  CBANK  ,  CFRNP,  CTHERM,  CTHPED, 

MAXMIN,  NUTRON,  RD,  RECAL,  SNGLES  ,  XTT 


Figure  4.  9.  Flowchart  for  subroutine  BANKR . 
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Subroutine  Name: 


Function: 


Calling  Program: 
Subroutines  Called: 
Labeled  COMMON: 


CHARGE  (Fig.  4.  10) 

Performs  photon  transport  calculations  in  air- 
over-ground  geometry.  Computes  and  scores 
energy  deposition  and  Compton  current  at  photon 
collision  sites.  Similar  to  TIGRE  program  used 
at  AFWL  for  prompt  gamma  EMP  drivers. 

BANKR,  GATR,  ROUTIN 

AZIRN,  GTISO,  KLEIN,  SEARCH 

CTHPED,  ENDXX,  MAXMIN,  RD,  RECAL,  XTT, 
ZSTOR 


S < 


Subroutine  Name: 


Function: 


Calling  Program: 
Subroutine  Called: 
Input  Parameters: 


Output  Parameters: 


Labeled  COMMON: 


DEP  (VOUTD,  WT,  PP,  TSD,  X2,  Y2,  Z2, 

T2,  VO)  (Fig.  4.  11) 

Computes  and  scores  energy  deposition  due  to 
thermal  neutron  collisions  and  nonthermal 
neutrons  treated  analytically. 

BANKR 

SEARCH 

VOUTD:  Velocity  of  neutron  after  energy  deposition. 

WT:  Weight  of  neutron. 

TSD:  Slowing  down  time  for  neutron 

X2,  Y2,  Z2:  Coordinates  of  the  neutron  energy 
deposition  position. 

T2:  Neutron  time  at  energy  deposition. 

VO:  Velocity  of  neutron  before  energy  deposition. 

WT:  Neutron  weight  after  energy  deposition. 

PP:  Gamma  weight. 

CFRNP,  CROSS,  CTHPED,  ENDXX,  MAXMIN, 

RD,  RECAL,  ZSTOR 


4. 11.  Flowchart  for  subroutine  DEP 


Subroutine  Name: 
Function: 

Calling  Program: 
Subroutines  Called: 
Labeled  COMMON: 


DOSEN  (Fig.  4. 12) 

Computes  and  scores  ionization  due  to  fast 
Monte  Carlo  neutron  collisions. 

NIES 

NINCS,  SEARCH,  XION 

CNXTAP,  ENDXX,  ENGNOS,  GROSS,  MAXMIN, 
PILL,  RD,  RECAL,  XTT,  ZSTOR 
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Figure  4.  12.  Flowchart  for  subroutine  DOSEN  . 


Subroutine  Name: 
Function: 

Calling  Program: 
Subroutines  Called: 
Input  Parameters : 


> 


| 

{ 


GATR  (NTIMES,  IFA)  (Fig.  4. 13) 

Determine  energy  of  capture  gammas.  Tracks 
analytic  thermal  gammas. 

BANKR,  XINIT 

CHARGE,  GTISO,  SEARCH 

NTIMES:  Number  of  gammas  to  be  produced. 
IFA  /  0:  Determine  only  gamma  energy. 


Labeled  COMMON:  CROSS,  C THERM,  CTHPED,  ENDXX,  MAXMIN, 

RD,  RECAL,  XTT,  ZSTOR 


; 


Figure  4. 13.  Flowchart  for  subroutine  GATR  . 
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Subroutine  Name: 
Function: 

Calling  Program: 
Labeled  COMMON: 


INELIN  (Fig.  4.14) 

Read  input  data  cards  related  to  inelastic 
scattering. 

INPUT 

INELC 
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Subroutine  Name: 

Function: 


Calling  Program: 
Subroutines  Called: 
Input  Parameters: 


Output  Parameters: 
Labeled  COMMON: 


KINNY  (SPDSQ,  U,  V,  W,  WATE,  NMED, 
LELEM)  (Fig.  4. 15) 

4 

Treats  the  inelastic  collisions,  determining 
final  speed  of  neutron  after  collision.  It  also 
determines  the  velocity  components,  U,  V, 
and  W. 

INELAS  which  is  called  by  SCATR. 

AB,  DCS,  ELAW,  LEVELS 

SPDSQ:  Neutron  speed  squared. 

U,  V,  W:  Velocity  components. 

WATE:  Weight  before  collision. 

NMED:  Medium  collision  took  place  in. 

LELEM:  Nuclide  collided  with. 

SPDSQ:  Neutron  speed  square  after  collision. 
V,  W:  Velocity  components. 

INELC 
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Figure  4. 15.  Flowchart  for  subroutine  KINNY. 

248 


4 


Subroutine  Name: 
Function: 

Calling  Program: 
Subroutines  Called: 
Input  Parameter: 


Labeled  COMMON: 


NIES  (NROUTE)  (Fig.  4.16) 

Zeros  scoring  arrays,  calls  XINIT,  ROUTEN, 
and  DOSEN. 

BANKR 

DOSEN,  ROUTIN,  XINIT 

NROUTE:  Specifies  what  action  is  to  be  taken 
by  NIES.  NROUTE  -  1  implies  start  of  program 
execution,  read  XINIT  input,  initialize  certain 
arrays  and  variables.  NROUTE  =  2  implies 
collision  information  is  to  be  analyzed.  NROUTE 
=  3  implies  end  of  an  energy  group. 

CALLE,  CBANK,  CNXTAP,  ENDXX,  MAXMIN, 
ONETWO,  PILL,  RD,  RECAL,  XTT,  XUNITS. 


Subroutine  Name: 
Function: 

Calling  Program: 
Subroutines  Called: 
Input  Parameters: 


Output  Parameters: 


Labeled  COMMON: 


RAGE  (AX,  AY,  Z,  T,  TAU,  Cl,  C2,  C3,  TS)  (Fig.  4. 17) 

Determine  new  epithermal  neutron  position  and  time, 
for  analytically  treated  neutrons. 

BANKR 

GTISO 

AX,  AY,  Z:  Old  neutron  position  coordinates. 

T:  Old  neutron  time. 

TS:  Slowing  down  time. 

TAU:  Fermi  age. 

AX,  AY,  Z:  New  neutron  position  coordinates. 

T:  New  neutron  time. 

Cl,  C2,  C3:  Components  of  an  isotropic  unit  vector. 


1 


CROSS 


i 

i 


;i 
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Subroutine  Name: 
Function: 

Calling  Program: 
Subroutines  Called: 
Input  Parameters: 

Output  Parameters: 

Labeled  COMMON: 


m1"  W-  \* '.v-', 1 1" *■.■*?* r 


RDIFF  (XO,  YO,  ZO,  TO,  NMED,  XI,  Yl,  Zl,  Tl)  (Fig. 4. 18] 
Determine  final  thermal  neutron  position  and  time. 

BANKR 
GTISO 

XO,  YO,  ZO:  Old  neutron  position 
TO:  Old  neutron  time 
NMED:  Medium  neutron  is  in 

XI,  Yl,  Zl:  Final  neutron  position 
Tl:  Final  neutron  time 

CROSS,  CTHPED 


i 

j 

i 

! 


i 
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Subroutine  Name: 
Function: 


Calling  Program: 
Subroutines  Called: 
Labeled  COMMON: 


ROUTIN  (Fig.  4.  19) 

Creates  photons  resulting  from  high-energy 
neutron  collisions,  giving  their  energies 
and  weights. 

NIES,  XINIT 

CHARGE,  NINCS 

CNXTAP,  CROUT,  GROSS,  RD, 

RECAL,  XTT 


255 
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Figure  4. 19.  Flowchart  for  subroutine  ROUTIN. 
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Subroutine  Name: 


Function: 

Calling  Program: 
Input  Parameters: 


Output  Parameter: 


SEARCH  (T,  TT,  N,  M,  K)  (Fig.  4.  20) 

SEARCH  is  used  to  determine  the  indexing 
values  which  specify  the  particular  array 
locations  in  which  information  is  to  be  stored. 

CHARGE,  DOSEN 

T:  A  value  within  the  range  TT(1)  to  TT(N). 

TT:  An  array  of  size  N  of  boundary  values 
in  increasing  numerical  order. 

N:  Size  of  TT  array. 

M:  Smallest  power  of  2  greater  than  or 

equal  to  N. 

K:  The  index  value. 
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Subroutine  Name:  SOURCE  (SPDSQ,  U,  V,  W,  X,  Y,  Z,  WATE, 

N,  NMEN,  NMED,  NREG)  (Fig.  4.  21) 

Function:  Specifies  the  initial  parameters  of  a  neutron. 

These  parameters  are:  neutron  energy,  direc¬ 
tion  of  its  first  flight,  and  its  statistical  weight. 
The  starting  coordinates  are  (0. ,  0. ,  0. ),  set 
by  a  data  card. 

Calling  Program:  MSOUR,  05RNIES 

Output  Parameters:  SPDSQ:  The  neutron  speed  squared. 

U,  V,  W:  The  direction  cosines  of  the 
source  neutron. 

WATE:  The  starting  weight  to  be  as¬ 

signed  to  the  source  neutron. 

Parameters  Not 

Set  or  Used:  X,  Y,  Z:  The  starting  coordinates  of  the 

source  neutron. 

N:  Integer  associated  with  a 

particular  neutron. 

NMEN:  Total  number  of  neutrons  to  be 

given  source  data. 

NMED:  The  medium  in  which  X,  Y,  and 

Z  lie. 

NREG:  The  region  in  which  X,  Y,  and 

Z  lie. 


Labeled  COMMON: 


CSOURC,  X UNITS 


Figure  4.  21.  Flowchart  for  subroutine  SOURCE. 


Subroutine  Name: 
Function: 

Calling  Program: 
Input  Parameter: 

Labeled  COMMON: 


SSTART  (NQST)  (Fig.  4.  22) 

Checks  to  see  that  correct  results  tape  is  mounted 
and  positions  the  desired  tape. 

05RNIES 

NQST:  Counter  used  to  indicate  how  many  energy 
bands  remain  to  be  processed  for  current  RUN. 

CSOURC,  CSTOR,  ENDXX,  RECAL,  XUNITS 
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Subroutine  Name : 

Function: 

Calling  Program: 

Subroutines  Called: 

Input  Parameters:  AZ:  Distance  above  or  below  air-ground  interface 

VIN:  Neutron  speed  at  start  of  Fermi  age  process 
NMED:  Medium  neutron  is  in 

Output  Parameters:  TAU:  Maximum  Fermi  age  decrement 

VOUT:  Neutron  speed  at  end  of  Fermi  age  process 

TS:  Slowing  down  time  corresponding  to  maximum 
Fermi  age  decrement 

TAUD:  Fermi  age  decrement  used  for  energy 
deposition  calculations 

VOUTD:  Neutron  speed  used  for  energy  deposition 
calculations 

TSD:  Slowing  down  time  used  for  energy  deposition 
calculations 

IFB:  =  1:  Cannot  treat  epithermal  neutron  analytically 
=  2:  Can  treat  epithermal  neutron  analytically 

Labeled  COMMON:  CTHPED 


TAUM  (AZ,  VIN,  NMED,  TAU,  VOUT,  TS,  TAUD 
TSD,  IFB)  (Fig.  4.  23) 

Determine  Fermi  ages,  slowing  down  times  and 
energy  for  analytically  treated  neutrons. 

BANKR 

SEARCH 
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Figure  4.  23.  Flowchart  for  subroutine  TAUM. 


Subroutine  Name:  UPDATE  (IFA)  (Fig.  4.  24) 

Function:  Writes  energy  band  data  and  other  information 

on  results  tape.  Determines  when  new  results 
tape  should  be  mounted.  Prints  out  counter 
arrays. 

Calling  Program:  05RNIES 

Input  Parameters:  IFA  =  1,  fast  results 

2,  thermal  results 

Labeled  COMMON:  COECS,  CSOURC,  CSTOR,  CVOL,  ENDXX, 

RECAL,  VTTX,  XUNITS,  ZSTOR 
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UPDATE 
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Subroutine  Name: 


XINIT  (NIST)  (Fig.  4.  25) 

Function:  Reads  and  prints  analysis  input.  It  also 

determines  constants  used  in  the  analysis 
routines  based  on  input  values. 

Calling  Program:  NIES 

Subroutine  Called:  TAUM,  GATR,  ROUTIN 

Output  Parameter:  NIST:  dummy 

Labeled  COMMON:  CBANK,  CNXTAP,  COECS,  CROUT, 

CTHERM,  CVOL,  ENDXX,  MAXMIN, 
ONETWO,  RECAL,  VTTX,  XUNITS, 
ZSTOR 


Read  problem  input  parameters 
for  analysis  subroutines. 


Calculate  bin  widths,  limits,  volumes  and 
other  analysis  subroutines  parameters. 
Also  modify  units  of  times  and 
radial  bin  boundaries. 


Figure  4.  25.  Flowchart  for  subroutine  XINIT. 


4.  6.  3.  2  Labeled  Common  Blocks 

In  this  subsection  we  define,  briefly,  all  of  the  variables  which 

are  to  be  found  in  labeled  common  blocks  of  the  analysis  section  of  OSRNIES. 

1.  COMMON/C  ALLE/ALLE 

ALLE  -  Energy  lost  due  to  neutrons  escaping  from  the  system 

2.  COMMON/ CBANK/ TIMMAX,  EPICUT,  ICOLL  (1000),  RREPI, 

RRTHM 

TIMMAX  -  Local  time  past  which  neutrons  no  longer  followed 

EPICUT  -  Epithermal  cutoff  in  (cm/sec)^ 

ICOLL  -  Array  containing  the  number  of  collisions  undergone 

by  the  given  neutrons 

RREPI  -  Russian  roulette  probability  to  determine  if  a  gamma 
ray  should  be  produced  from  an  epithermal  neutron 
collision 

RRTHM  -  Russian  roulette  probability  to  determine  if  a  gamma 
ray  should  be  produced  from  a  thermal  neutron 
collision 

3.  COMMON/CFRNP/FRNP 

FRNP  -  Fraction  of  the  absorption  cross  section  for  air 

which  results  in  photon  emission  for  epithermal 
and  thermal  neutron  collisions 

4.  COMMON/CNXTAP/NXTAP,  SIG  (1700),  NPQ,  ISTRT  (40), 

EBOTI,  ETOPI,  NPT  (40) 

NXTAP  -  Gamma  production  cross  section,  Q,  and  energy 

deposition,  data  file  number 

SIG  -  Array  into  which  gamma  production  cross  section 

Q,  and  energy  deposition,  data  are  stored 


NPQ 

ISTRT 


Record  length  of  the  above  (SIG)  data 


Array  containing  the  starting  locations  of  the 
individual  cross  section  in  the  SIG  array 

2 

EBOTI  -  Lower  energy  (cm/sec)  boundary  of  the  SIG  data 

2 

ETOPI  -  Upper  energy  (cm/sec)  boundary  of  the  SIG  data 

NPT  -  Array  containing  the  size  of  the  individual  cross 

sections  in  the  SIG  array 

COMMON/COECS/I1,  12,  13,  N2 

11  -  Maximum  number  of  spatial  bins 

12  -  Maximum  number  of  time-spatial  bins 

13  -  2*  (11  +  12) 

.N2  -  Number  of  items  to  be  stored.  Possible  items  are: 

ionization  due  to  neutrons,  ionization  due  to  gam¬ 
mas,  horizontal  (or  radial)  currents,  and  vertical 
(or  theta)  currents. 

COMMON/CROSS/GSIGA,  GD,  RTPI,  PI,  EPS,  ASIGA,  ADXX 
Thermal  absorption  cross  section  in  ground 
Thermal  diffusion  coefficient  in  ground 
-  yfTT 

7T 

e  =  2.  718282 

Thermal  absorption  cross  section  in  air 


GSIGA 

GD 

RTPI 

PI 

EPS 

ASIGA 


ADXX 


Thermal  diffusion  coefficient  in  air 


7.  COMMON/CROUT/NPHOT,  WPHOT,  WTLIM,  TENWTL 

NPHOT  -  Number  of  photons  produced  by  each  inelastic 
neutron  collision 

WPHOT  -  Floating  point  representation  of  NPHOT 

WTLIM  -  Weight  below  which  Russian  Roulette  is  played 
in  ROUTIN  to  determine  if  photon  is  followed 

TENWTL  -  10  *  WTLIM,  used  in  determining  survival 

probability  of  photon  in  ROUTIN 

8.  COMMON/ CSOURC/IFAG,  INBA,  IRESLT,  NBA,  NEB, 

NEBL,  NRPT,  NOREEL,  NQUIS 

IF  AG  -  =1,  start  new  results  tape 

INBA  -  Energy  band  currently  processing 

IRESLT  -  Unit  number  of  results  tape 

NBA  -  Number  of  runs  (batches) 

NEB  -  Number  of  energy  bands 

NEBL  -  Number  of  energy  bands  left  to  be  processed 

for  current  runs 

NRPT  -  Number  of  runs  per  results  tape 

NOREEL  -  Reel  number  of  current  results  tape 
NQUIS  -  Number  of  runs  left  to  be  processed 
E  -  Energy  band  boundaries  (monotonically  decreasing) 

9.  COMMON/CSTOR/ID 

ID  -  Run  ID  number 
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10.  COMMON/C THERM/NTHERM,  TMULT(2),  A3 (2) 


NTHERM  -  Number  of  gammas  produced  by  thermal  neutrons 
on  the  terminating  collisions 

TMULT  -  Neutron  capture  gamma  multiplicities:  TMULT(l)- 
air,  TMULT(2)-ground 

A3  -  Boundary  value  (cm)  to  determine  whether  Monte 

Carlo  or  analytic  method  should  be  used  to  track 
thermal  neutrons  A3(l)-air  (elevation),  A3  (2)- 
ground  (depth) 

11.  COMMON/C THPED/THPE ED,  THPSQ 

THPEED  -  Velocity  of  thermal  neutrons  (cm/sec) 

THPSQ  -  THPEED2  (cm/sec)2 

12.  COMMON/C VOL/DELT  (41),  VOLB  (360) 

DELT  -  Time  bin  widths  (sec) 

3 

VOLB  -  Volumes  of  the  spatial  bins  (m  ) 

13.  COMMON/ ENDXX/ENDEP  (10400),  GENDEP  (10400), 

XCR  (10400),  XCT  (10400),  IMAX,  JMAX,  KMAX 

ENDEP  -  Array  in  which  energy  deposition  due  directly  to 
neutrons  is  stored 

GENDEP  -  Array  in  which  energy  deposition  due  to  gammas 
created  by  neutron  collisions  is  stored 

XCR  -  Array  in  which  radial  or  horizontal  current 

component  is  stored 

XCT  -  Array  in  which  theta  or  vertical  current  component 

is  stored 
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14.  COMMON/ ENGNOS/WO(28),  WN(26) 

WO  -  Energy  levels  (Q's)  of  the  discrete  inelastic 

oxygen  scatterers 

WN  -  Energy  levels  (Q's)  of  the  discrete  inelastic 

nitrogen  scatterers 

15.  COMMON/ GROSS/NDISC (2),  NCONT(2),  GAMENG(2, 30) 

NDISC(I)  -  Number  of  discrete  photon  energy  cross  sections 
for  medium  I 

NCONT(I)  -  Number  of  continuum  photon  energy  cross 
sections  for  medium  I 

GAMENG 

(I,  J)  -  Energy  (MeV)  of  photon  resulting  from  cross 

section  J  in  medium  I 


16.  COMMON/INELC/SP,  UP,  VP,  WP,  E,  S,  NINLS,  IN,  NR(120), 

Q(l,  120),  NE(1, 120),  EL  (1,1,1),  SIGL  (1,1,1),  NU(120),  ELAW 
(1,120),  NLAW  (1,120),  NENG(120),  ENG(5,2),  NETAB1(5,2) 

ETAB1(23,  5, 2),  PTAB1(23,  5, 2),  NETAB2(1),  ETAB2 (1,1), PTAB2 (1,1), 
FM(1 , 120),  EA(1),  EB(1),  TA(120),  TB(120),  TC(120),  NAPP 
(120),  NMEDE(2, 120),  NELE(2,120),  NAE(120) 


SP 

UP 

VP 

WP 

E 

S 

NINLS 


Velocity  of  resultant  neutron  (cm/sec) 
X  axis  unit  component  of  SP 
Y  axis  unit  component  of  SP 
Z  axis  unit  component  of  SP 
Energy  of  incident  neutron  (eV) 
Velocity  of  incident  neutron  (cm/sec) 
Number  of  inelastic  scatterers 


IN 


Scatterer  involved  in  neutron  collision 


16.  (Cont.) 
NR 


The 


Q 

NE 

E 

SIGL 

NU 

ELAW 

NLAW 


scatter er  flag,  =1  for  discrete  levels 

Q  value  for  scatterer 

dummy 

dummy 

dummy 

Scatterer  flag,  «1  non-discrete  scatterer 

Energy  threshold  for  non-discrete  scatterers 

Indicates  distribution  to  be  used  for  selecting 
i  esultant  neutron  energy  6 


following  arrays  deal  with  tabulated  E*  vs  P(E*)  scatterers 
NENG  -  Number  of  energy  levels 

ENG  -  Energy  levels  of  incident  neutron 


NETAB1 


Number  of  individual  P(E*)  values  for  given 
incident  energy 


ETAB1 

PTAB1 


Resultant  neutron  energy  associated  with  a 
particular  cummulative  P(E*) 

Cummulative  P(E*)’s 


NETAB2  -  dummy 
ETAB2  -  dummy 
PTAB2  -  dummy 
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16.  (Cont. ) 

FM  -  Multiplicity  associated  with  particular  non 

discrete  scatterer 

EA  -  dummy 

EB  -  dummy 

TA  -  Maxwellian  coefficient 

TB  -  Maxwellian  coefficient 

TC  -  Maxwellian  coefficient 

NAPP  -  Number  of  media  scatterer  is  in 

NMEDE  -  Medium  number  scatterer  is  in 

NELE  -  Scatterer  number  with  respect  to  media 

NAE  -  Tabulated  E*  vs  P(E*)  scatterer  number 

17.  COMMON/M AXMIN/RMAX,  RMIN,  TMAX,  ZMAX,  ZMIN, 
RMAXZ,  A 

RMAX  -  Maximum  radius  of  interest  (cm) 

RMIN  -  Minimum  radius  of  interest  (cm) 

TMAX  -  Maximum  time  of  interest  (sec) 

ZMAX  -  Maximum  Z  value  of  interest  (cm), 

source  at  Z  =  0. 0. 

ZMIN  -  Minimum  Z  value  of  interest  (cm) 

RMAZ  -  RMAX  *  RMAX 

Height  of  burst  of  source  (cm) 


A 


18.  COMMON/ONE  TWO/IRP,  ISP,  IWANT 


IRP  -  dummy 

ISP  -  dummy 

IWANT  -  Flag  which  specifies  what  collisions  are  to 
be  analyzed. 

19.  COMMON/PILL/XWT,  EN,  EA 

XWT  -  Weight  of  unbiased  neutron 

EN  -  Neutron  energy  (MeV)  before  collision 

EA  -  Neutron  energy  (MeV)  after  collision 

20.  COMMON/RD/NCOLL,  NAME,  SPDSQ,  U,  V,  W,  X,  Y,  Z,  WATE, 
SPOLD,  UOLD,  VOLD,  WOLD,  XOLD,  YOLD,  ZOLD,  OLDWT, 
THETM,  PSTL,  ETAUSD,  NGROUP,  LELEM,  NR  EG.  NMED. 
NAMEX,  WATEF,  BLZNT,  BLZ,  LAMBDA,  SO,  SI,  ETATH, 
FONE,  EXTRA1,  EXTRA2 

The  05R  collision  parameters  associated  with  each  neutron 
collision 

21.  COMMON/RECAL/PCOUNT(30),  ICOUNT(30),  ZCOUNT(30), 
KCOUNT(30),  IPKT,  DKT,  EPICUT 

PCOUNT  -  Floating  point  counters  used  to  check  energy 
conservation,  energy  of  neutrons  and  gammas 
produced  in  a  given  energy  group 

ICOUNT  -  Integer  counters  used  to  indicate  number  of 
neutrons  and  gammas  produced,  and  how 
many  times  certain  types  of  collisions  occur 
in  a  given  energy  group 

ZCOUNT  -  Accummulative  PCOUNT  over  a  given  run 
Accummulative  ICOUNT  over  a  given  run 


KCOUNT 


21. 

(Cont  ) 

IPKT 

Size  of  PCOUNT  and  ZCOUNT  arrays 

IKT 

Size  of  ICOUNT  and  KCOUNT  arrays 

EPICUT 

Epithermal  cutoff 

22. 

COMMON/VTTX/NBATCH,  NHLS 

NBATCH 

Number  of  batches  per  energy  group 

NHIS 

Number  of  neutron  histories  pe/'  batch 

23. 

COMMON/XTT/T(l 000),  PP,  EG 

T 

Elapsed  time  of  neutrons 

PP 

Weight  of  phoion 

EG 

Energy  of  photon 

24. 

COMMON/XUNITS/IP,  JP 

IP 

Card  input  unit  number 

JP 

Printed  output  unit  number 

25. 

COMMON/ZSTOR/RADC(41),  NRA,  MRAY,  TIMC(41),  NTIM 
MTYM,  ZADC  (41),  NZA,  MZAY,  RADCS(41),  ANGC(21), 
NAN,  MAN,  NCOSA,  NR  ADA,  NTIMA,  NZMA 

RADC 

Radial  bin  boundaries  (cm) 

NRA 

NR  ADA  +  1 

MRAY 

-  MOWER  (NRA) 

TIMC 

Time  bin  boundaries  (sec) 

NTIM 

-  NTIMA  +  1 
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(Cont. ) 

MTYM 

-  MOWER  (NTIM) 

ZADC 

Horizontal  bin  boundaries  (cm) 

NZA 

-  NZMA  +  1 

MZAY 

-  MOWER  (NZA) 

RADCS 

-  RADCS  (I)  =  RADC  (I)  *  RADC  (I) 

angc 

Angular  bin  boundaries  (cos  6) 

NAN 

-  NCOSA  +  1 

MAN 

-  MOWER  (NAN) 

NCOS  A 

Number  of  angular  bins 

NR  ADA 

Number  of  radial  bins 

NTIMA 

Number  of  time  bins 

NZMA 

Number  of  horizontal  bins 
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4.  6.  3.  3  External  Tapes  Required 

As  the  program  is  currently  run  three  external  tapes  are 
used,  with  only  one  being  required  at  any  one  time.  The  program  file  tape 
is  initially  requested  and  recreated  on  an  internal  file  being  modified  as 
desired.  Next,  the  data  tape  containing  files  NPTAPE  (file  4),  NCONTP 
(file  3),  NXTAP  (file  14),  and  TAPE  15  (file  15)  is  mounted.  After 
copying  these  files  to  internal  files  it  may  be  removed.  Finally,  the 
results  file  ERESLT  (TAPE  16)  is  mounted.  As  required,  extra  reels  for 
TAPE  16  are  mounted. 

4.  6.  3.  4  Data  Files  Used 

NPTAPE:  PHI  tape.  Generated  by  a  COD678  CODE  8  run. 
Contains  data  necessary  to  employ  the  Coveyou  technique  for  anisotro,  ■. 
angular  scattering.  See  Ref.  ( 12)  for  further  details. 

NCONTP:  05R  system  data  tape.  Generated  by  a  CCu67ii 
CODE  6  run.  Contains  cross  section  data  for  05R  neutron  collision 
routines.  See  Ref.  ( 12)  for  further  details. 

NXTAP:  The  file  NXTAP  contains  nuclear  data  which  are  re¬ 
quired  for  the  analysis  of  the  neutron  histories,  but  which  are  not  available 
on  the  05R  system  tapes.  The  data  are  read  as  a  separate  record  for 
each  energy  supergroup.  For  each  energy  supergroup,  the  complete  array 
contains  several  types  of  nuclear  data.  First  is  cross  section  data  for 
gamma  production  by  neutron  inelastic  scattering  in  air.  These  data  are 
in  the  form  of  cumulative  probability  distributions.  Next  is  similar  data 
for  ground.  The  gamma  production  data  are  for  a  particular  ground  com¬ 
position,  and  must  be  changed  whenever  a  different  composition  is  to  be 
investigated.  These  data  are  followed  by  a  set  of  neutron  kerma  values  for 
ground.  These  are  used  to  compute  the  energy  deposition  in  ground. 

This  is  followed  by  a  set  of  Q  values,  as  discussed  in  Section  VI.  5  of 


279 


Ref.  (  1  ).  The  values  on  tape  are  normalized  by  the  total  cross  section; 
i.  e. , 

<3 - LJ_ 

so  that,  the  quantity  (W^  -  Wa)  (E  +  Q)  is  scored  at  every  collision  to  account 
for  energy  deposition  by  neutron  absorption.  Finally  there  follows  a  set 
of  values  for  the  average  energy  transfer  available  to  electrons  from 
elastic  scattering  neutron  collisions  in  air.  These  are  the  data  shown  in 
Fig.  3  of  Ref.  (1  ),  and  are  used  to  score  ionization  in  air. 

TAPE  15:  The  tape  15  was  created  to  eliminate  a  large  number 
of  data  cards  from  the  input  decks,  these  data  being  the  same  for  all  runs. 
The  binary  file  consists  of  six  records,  five  of  which  are  used  in  TAUM 
the  other  in  ROUT  IN. 

RECORD  1:  NVA,  NVG 

NVA  -  Number  of  energy  groups  represented  by 
data  in  TAUS  in  air. 

NVG  -  Number  of  energy  groups  represented  by 
data  in  TAUS  in  ground. 

RECORD  2:  (VSTEP  (I,  1)  I  =  1,  NVAP1) 

NVAPI  =  NVA  +  1 

NSTEP:  Energy  group  boundaries  (v^),  in 

increasing  value  (vi  <  v\  +  ^),  in  cm/ sec. 

RECORD  3:  (  TAUS  (I,  J,  1),  J  =  1,  2),  I  =  1,  NTAUSA) 

J  =  1,  Fermi  ages  (f.  .)  for  air 

J 

J  =  2,  Slowing  down  time  (t^  for  air 
NTAUSA  =  (NVAPI  *  NVA)/2 
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f .  .  -  is  the  Fermi  age  for  going  from  v,  to  v . 

J  K  J 

where  v,  >  v  .  . 
k  >  J 

t.  .  -  is  the  slowing  down  time  corresponding  to 

k,  j 

fk,  j  • 

Then,  the  f’s  are  stored  in  the  following  manner: 
f  •  f  •  f  -f  ■  f  •  f 

l2 ,  1  ’  1 3,  2  ’  1 3,  1  ’  14,  3  ’  *4,  2  »  *4,  1  ’ 

;  f NVA,  NVA  -  1  5  ’  f NVA, 1  5 

f  NVA  +  1,  NVA  ;  ;  fNVA+  1,  1. 

The  tfs  are  stored  in  a  similar  manner. 

RECORD  4:  (VSTEP  (I,  2)  I  =  1,  NVGP1  ) 

Same  as  record  2  except  it  is  for  the  ground 

RECORD  5:  (TAUS  (I,  J,  2),  J  =  1,  2  ),  1=1,  NTAUSG) 

Same  as  record  3  except  it  is  for  the  ground. 


RECORD  6:  NSENG,  NEENG,  ENG,  PGE,  GE 

DIMENSION  NSENG  (3),  NEENG  (3),  ENG  (21), 

PGE  (45,  21),  GE  (45,  21) 

Three  gamma  production  cross  sections  in  file  14  have  the 
gamma  energy  in  the  form  of  tabulated  probability  distributions  which  are 
represented  in  the  above  arrays. 

NSENG  (I)  -  Starting  energy  index  value  for  the  Ith  cross  section. 

th 

NEENG  (I)  -  Final  energy  index  value  for  the  I  cross  section. 

ENG  (J)  -  Upper  energy  value  (MeV)  for  the  incident  neutron 
th 

which  belongs  to  the  J  probability  distribution. 

th 

(PGE  (K,  J),  K  =  1,  45)  -  J1  accumulative  probability  distribution. 
(GE  (K,  J),  K  =  1,  45)  -  corresponding  gamma  energy  values. 


IRESLT:  05RNIES  results  file.  Each  reel  of  IRESLT 
consists  of  an  identification  record  followed  by  sets  of  data  records, 
each  set  consisting  of  eight  records.  These  are  NRPT  *  NEB  data  sets, 
where  NRPT  is  the  number  of  batches  per  results  tape  reel  and  NEB  is 
the  number  of  energy  groups  per  batch.  The  first  record  consists  of: 

ID:  Run  identification  number. 

NOREEL:  IRESLT  tape  reel  number. 

NZ:  Number  of  items  stored. 

NEB:  Number  of  energy  groups . 

NBA:  Number  of  batches . 

NTIMA:  Number  of  time  bins. 

NRADA:  Number  of  radial  bins. 

NCOSA:  Number  of  angular  bins. 

NZMA:  Number  of  horizontal  bins. 

TIMC:  Time  bin  boundaries. 

RADC:  Radial  bin  boundaries. 

ANGC:  Angular  bin  boundaries. 

ZADC:  Horizontal  bin  boundaries. 

E:  Energy  bin  boundaries. 

NRPT:  Number  of  batches  per  reel  for  IRESLT  file. 

NHIS:  Number  of  neutrons  per  energy  band  run. 

DELT:  Time  bin  widths. 

VOLB:  Volume  of  the  spatial  bins. 

Each  record  in  a  data  set  consists  of  one  of  the  arrays 
dimensioned  10400  which  are  found  in  the  common  block  ENDXX.  The 
results  are  stored  in  these  arrays  as  if  they  were  dimensioned  (NTIMA, 
NRADA,  NCOSA  +  NZMA).  The  eight  records  are: 
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ENDEP 

GENDEP 

XCR 

XCT 

ENDEP 

GENDEP 

XCR 


Non-thermal 


Thermal 


4.  7  INSTRUCTIONS  FOR  USING  THE  05NPT  CODE 

4.  7. 1  Q5NPT  Input  Instructions 

The  following  is  a  description  of  the  input  required  for  05NPT. 

CARD  A:  TITLE 

FORMAT  (10A6) 

TITLE  (dimensioned  10):  A  60  character  title  which  will 
appear  on  all  plots. 


CARD  B:  MRUN,  IRTT,  N2,  IFGTP,  IFGTIO,  IFGTIP,  IFGINO,  IFGINP 
FORMAT  (1215) 

MRUN  -  Number  of  batches  data  represent 

IRTT  -  Data  set  on  file  KP  which  contains  the  results 
for  MRUN  batches 


N2  -  =  4,  which  is  the  number  of  items  written  on  the 

energy  band  results  tape 

IFGTP  (dimensioned  5): 

IFGTP(l)  =  0,  plot  and/or  print  out  neutron  plus  gamma 
ionization  results 

IFGTP (2)  =  0,  plot  and/or  print  out  neutron  ionization 
results 

IFGTP(3)  =  0,  plot  and/or  print  out  gamma  ionization 
results 


IFGTP(4)  =  0,  plot  and/or  print  out  radial  and  horizontal 
currents 

IFGTP(5)  =  0,  plot  and/or  print  out  theta  and  vertical 
currents 


IFGTIO 

IFGTIP 

IFGINO 

IFGINP: 


Time-dependent  results  are  to  be  printed  out 
Time-dependent  results  are  to  be  plotted 
Time- independent  results  are  to  be  printed  out 
Time- independent  results  are  to  be  plotted 
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CARD  C: 


CARD  D: 


CARD  E: 


CARD  F’s: 


SCALE 

FORMAT  (5E10.  3) 

SCALE  (dimensioned  5):  Scale  factors  by  which  the  different 
items  2 re  to  be  multiplied. 

NIA(I),  N3A(I),  N4A(I),  A3A(I),  A4A(I),  DXA(I),  XMAXA(I), 
XMINA(I) 

FORMAT  (315,  5X,  5E10.  3) 

5  cards,  1  for  each  item  on  results  tape 
N1A(I)  -  Graph  type  for  item  I 

1  -  X-axis  log,  Y-axis  log 

2  -  X-axis  lin,  Y-axis  log 

3  -  X-axis  lin,  Y-axis  lin 

4  -  X-axis  log,  Y-axis  lin 

N3A(I)  -  Number  of  X-axis  decades  if  logarithmic 

N4A(I)  -  Number  of  Y-axis  decades  if  logarithmic 

DXA(I)  -  X-axis  interval  between  tick  marks  on  graph 
if  logarithmic 

XMAXA(I):  Maximum  X  value  allowed  for  plotting 
XMINA(I):  Minimum  X  value  allowed  for  plotting 

IR1,  IRTT,  NR,  ISP  (read  in  routine  SUMOBD) 

FORMAT  (413) 

IR1:  Batch  number  to  start  processing  energy  band  result 
tape  with 

IRTT:  Data  set  on  file  KP  which  contains  the  results  for  IR1-1 
batches 

NR:  Batch  number  to  terminate  processing  with 

ISP:  =  1,  after  processing  each  energy  band  results  tape  a 
data  set  is  written  on  file  KP 

(P(I),  1=1,  NEB)  (read  in  routine  SUMOBD) 

FORMAT  (8F10.  3) 

P:  Probabilities  to  be  applied  to  the  energy  band  results 

so  as  to  simulate  a  given  spectrum 
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4.  7.  2  05NPT  Output 


File  6  is  used  for  output. 

First  page:  Summary  of  card  input 

Second  page:  VCOUNTS:  Integer  counters  normalized  to  per  neutron 
XCOUNTS:  Floating  point  counters  normalized  to  per 
neutron 

The  VCOUNT  and  XCOUNT  arrays  are  analagous  to  the 
KCOUNT  and  ZCOUNT  arrays  of  Section  4.  6.  2. 


All  subsequent  pages:  Normalized  results  in  tabular  form  containing 
value,  its  standard  deviation,  and  its  percent  error 


4.  7.  3  Structure  of  the  Q5NPT  Code 
4.  7.  3. 1  Flowcharts 

On  the  following  pages  are  descriptions  of  the  various  routines 
in  05NPT,  excluding  the  PLOT  package  routines.  Flowcharts  are  also 
presented  for  the  more  complex  routines.  Some  of  the  subroutines  called 
by  05NPT  are  shown  in  Fig.  4.  26. 


i 


t 


Main  Routine:  05NPT  (Figs.  4.  27-4.  28) 

Function:  Produces  the  desired  information  in  an  orderly 

fashion  calling  necessary  subroutines. 

Subroutines  Called:  OUTPUT,  PLOTER,  SUMOBD,  XASA 

Labeled  COMMON:  CKT,  SINGLS 


Read  in  input  specifying  number  of  runs,  plots  and  output  desired, 
and  graph  parameters. 


CALL 

SUMOBD 


Position  spectrum  results  tape  and  read  desired  data  onto  disk 
file. 

Set  up  X-axis  coordinates  for  plotting  data. 

Set  up  some  normalizing  variables. 

Normalize  counters. 


4) 

5) 


6) 


Set  some  spatial  bin  dependent  variables. 

Print  output  heading  if  printed  outjjut  desired. 

Determine  normalized  time  bin  values,  associated  standard 
deviations  and  percent  errors  for  all  time  bins  associated  with 
a  given  spatial  bin. 

Printed  output  if  time  dependent  output  for  spatial  bins  desired. 
Determine  normalized  spatial  bin  value,  associated  standard 
deviation,  and  percent  error. 

Print  above  item  5  values  if  time  dependent  output  desired. 


Figure  4.  27.  Flowchart  for  the  05NPT  program. 
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Figure  4.  28.  Flowchart  for  the  05NPT  Program. 
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Subroutine  Name: 


OUTPUT  (NI,  NZ,  INDX,  WORD,  VA,  SDA, 

VB,  SDB,  VC,  SDC,  X,  TITLE) 

Function:  Put  data  in  tabular  form  in  the  output. 

Called  By:  05NPT 

Input  Parameters:  Nl  -  X  -  axis  type  (e.  g. ,  local  time,  range,  etc. ) 

NZ  -  Number  of  points 

INDX  -  Flag  specifying  whether  data  is  result 
of  neutrons  and  gammas,  neutrons,  or 
gammas 

WORD  -  Specifies  type  of  current  (e.  g. ,  radial, 
theta,  etc. ) 

VA  -  Fast  data  array 

SDA  -  Standard  deviation  of  VA  array  data 

VB  -  Thermal  data  array 

SDB  -  Standard  deviation  of  VB  array  data 

VC  -  Fast  and  Thermal  data  array 

SDC  -  Standard  deviation  of  VC  array  data 

X  -  Independent  variable  data 

TITLE  -  Title  for  output  table 
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Subroutine  Name: 

Function: 

Called  By: 
Subroutines  Called: 
Input  Parameters: 


w*rr[*— nit^iT-T^r— 


PLOTER  (NT,  N2,  XA,  SA,  INDX,  IFY,  ITIME) 
(Fig.  4.29) 

Plot  desired  data  in  histogram  form  with  error 
bars  on  the  total  curve. 

05NPT 


GRAPH  and  other  plotting  routines 


NT 

N2 

XA 

SA 


X-axis  type  (e.  g.  local  time,  range,  etc. ) 
Number  of  points 
Histogram  X-axis  coordinates 
Error  bar  X-axis  coordinates 


INDX  -  Flag  specifying  whether  data  are 

result  of  neutrons  and  gammas,  neutrons, 
or  gammas. 

IFY  -  Flag  specifying  Y-axis  label 

ITIME  -  =1,  Time  dependent  data 
=2,  Time  independent  data 


Labeled  COMMON: 


CK1F  CMXMN,  CPLOT 


START 


1)  Set  some  subroutine  variables 

2)  Determine  Y-axis  coordinate  limits 
for  plot 

3)  Plot  total  curve 

4)  Draw  error  bars 

5)  Plot  fast  curve 

6)  Plot  thermal  curve 

7)  Label  gra^h 


RETURN 


Figure  4.  29.  Flowchart  for  subroutine  PLOTER 


3  mini  l  mm&mtm 


Subroutine  Name:  SUMOBD  (Fig.  4.  30) 

Function.  Create  file  KP  (ionization  rates  and  currents 

for  a  particular  neutron  spectrum)  using  results 
on  file  IRESLT  (ionization  rates  and  currents 
for  energy  bands). 

Called  By:  05NPT 

Subroutines  Called:  SUMS 

Labeled  COMMON:  CKT,  SINGLS,  C UNITS 
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Figure  4.  30. 


Flowchart  for  the  subroutine  SUMOBD 
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Subroutine  Name: 


SUMS  (IMAX,  JMAX,  KMAX,  SEMP,  ENDEP, 
ICOC) 

Function:  Creates  the  sum  over  the  batches  and  the  sum 

of  the  squares  over  the  batches 

Called  By:  SUMOBD 

Input  Parameter:  IMAX :  Number  of  time  bins 

JMAX:  Number  of  radial  bins 

KMAX:  Number  of  angular  plus  horizontal  bins 

ENDEP:  Array  of  values  to  be  added  to 

accumulative  total  over  the  batches 
ICOC:  Dummy 

Output  Parameter:  SEMP:  Array  containing  the  time  integrated 

values  of  the  ENDEP  array 

Labeled  COMMON:  SINGLS 
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Subroutine  Name: 
Function: 

Called  By: 

Input  Parameters: 


Output  Parameter: 


XASA  (AA,  AM,  IP,  NO,  XA,  SA,  KF) 

Determine  the  values  to  be  used  for  the  X-axis 
when  plotting  the  data. 

05NPT 

AA:  Bin  boundaries 

AM:  X-axis  coordinate  of  the  first  point  to 

be  plotted. 

IP:  =1  or  4  log  X-axis 

=2  or  3  lin  X-axis 

NO:  Number  of  bins 

KF:  =1,  increment  through  AA  array  by  1 

=  -1,  increment  through  A  A  array  by  -1 

XA:  Histogram  X-axis  coordinates 
SA:  Error  bar  X-axis  coordinates 
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4.  7.  3.  2  Common  Blocks 


In  this  sub-section  we  define,  briefly,  all  of  the  variables  which 
are  to  be  found  in  the  common  blocks  of  05NPT  and  subroutines  PLOTER, 
OUTPUT,  SUMOBD,  and  SUM. 

I.  COMMON/CKT/IIKT,  IPKT 

IIKT  -  Size  of  KCOUNT  and  ICOUNT  arrays 
IPKT  -  Size  of  PCOUNT  and  ZCOUNT  arrays 

n.  COMMON/  CKIF/KIF 

KIF  -  Indicates  the  type  of  line  used  in  plotting  a  set 
of  data 

KIF  -  0  - - - 

KIF  =  1  . . . . 

KIF  =  6 - 

n.  COMMON/ CMXMN/XMXL,  XMNL,  YMAX,  YMIN 

XMXL  -  Plots  X-axis  maximum 
XMNL  -  Plots  X-axis  minimum 
YMAX  -  Plots  Y-axis  maximum 
YMIN  -  Plots  Y-axis  minimum 


IV. 


C  OMMON/C  P  LOT/ YAP  (40), 
TITLE (10),  N1A(5),  N3A(5), 
XMAXA(5),  XMINA(5) 


YBP(40),  YCP(40),  SDCP(40) 
N4A(5),  A3A(5),  A4A(5),  DXA(5), 


YAP 

Y  coordinates  of 

YBP 

Y  coordinates  of 

YCP 

Y  coordinates  of 

plotted 

the  fast  data  to  be  plotted 
the  thermal  data  to  be  plotted 
the  fast  and  thermal  data  to  be 


SDCP 


Standard  deviation  of  the  fast  and  thermal  data 
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TITLE  -  A  60  character  title  to  be  printed  on  each  plot 

NIA  -  Specifies  type  of  X-axis  and  Y-axis  desired 

N3A  -  Specifies  number  of  decades  for  X-axis  if 

logarithmic 

N4A  -  Specifies  number  of  decades  for  Y-axis  if 
logarithmic 

A3A  -  Specifies  X-axis  plot  length 

A4A  -  Specifies  Y-axis  plot  length 

DXA  -  Specifies  X-axis  tick  mark  interval  for  linear 
X-axis 

XMAXA  -  Maximum  X  value  to  be  plotted 

XMINA  -  Minimum  X  value  to  be  plotted 

COMMON/C  UNITS/IP,  JP 

IP  -  Input  unit 

JP  -  Output  unit 

COMMON/SINGLS/KP,  IER,  IEW,  IERA,  IEWA,  IET,  IRESLT, 

II,  12,  13,  N2 

KP  -  Spectrum  results  file 

IER  -  Temporary  file  to  be  read  from,  contains  spec¬ 

trum  results 

IEW  -  Temporary  file  to  be  written  on,  contains  the 
latest  spectrum  results.  The  values  for  IER 
and  IEW  are  alternated. 

IERA  -  Temporary  file  to  be  read  from,  contains  the 
current  spectrum  run  sum 

IEWA  -  Temporary  file  to  be  written  on,  contains  the 
new  spectrum  run  sum.  The  values  for  IERA 
and  IEWA  are  alternated 

IET  -  Temporary  file  to  be  written  and  then  read. 

Contains  fast  results  which  are  to  be  added  to 
thermal  results  to  arrive  at  total  results 

IRESLT  -  Energy  band  results  file 

II  -  Size  of  the  F  and  G  arrays 
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12  -  Size  of  the  D  and  E  arrays 

13  -  2*  (II  +  12) 

N2  -  Number  of  items  of  results  data  stored  (i.  e. , 

(1)  Ionization  dose  due  to  neutron  ionizations; 

(2)  Ionization  dose  due  to  gammas  resulting  from 
neutron  collisions;  (3)  Radial  or  horizontal  cur¬ 
rents;  and  (4)  Theta  or  vertical  currents) 

VII.  COMMON  BSC (360),  BSSC(360),  ASC (10400),  ASSC (10400), 

ASA(10400),  ASSA(10400),  ASB(10400),  ASSB(10400),  DUM(10400), 
BSA(360),  BSSA(360),  BSB(360),  BSSB(360) 

In  the  main  routine  these  arrays  are  used  to  store  spectrum  re¬ 
sults  which  are  to  be  plotted  and/or  printed  out  after  slight  modification. 

In  the  subroutine  SUMOBD  the  first  4  arrays  are  used  to  store  information 
which  is  to  be  written  on  the  spectrum  results  file,  while  the  next  five  are 
used  to  store  information  read  from  the  energy  band  results  file.  The 
arrays  dimensioned  360  are  used  to  store  time  independent  data,  while 
those  dimensioned  10400  are  used  to  store  time  dependent  data. 


4.  7.  3.  3  05NPT  Data  File  Structure 


There  are  two  main  data  files  used  in  05NPT;  they  are  IRESLT 
and  KP.  IRESLT  is  the  input  data  file  and  KP  is  the  output  data  file,  along 
with  being  the  file  used  for  obtaining  the  information  to  be  plotted.  IRESLT 
has  been  fully  described  in  the  section  on  05RNIES  data  files.  The  file  KP 
consists  of  data  sett  which  contain  16  records  each.  Each  data  set  represents 
the  spectrum  results  for  a  given  number  of  batches.  The  first  record  con¬ 
sists  of:  ZCOUNT,  KCOUNT,  NR,  RADC,  TIMC,  ZADC,  ANGC,  NR  ADA, 
NTIMA,  NZMA,  NCOSA,  NBATCH,  NHIS,  DELT,  VOLB. 

ZCOUNT  (dimensioned  30):  Floating  point  counters 

KCOUNT  (dimensioned  30):  Integer  counters 

NR:  Number  of  batches  data  represents 

RADC:  Radial  bin  boundaries 

TIMC:  Time  bin  boundaries 

ZADC:  Horizontal  bin  boundaries 

ANGC:  Angular  bin  boundaries 

NRADA:  Number  of  radial  bins 

NTIMA:  Number  of  time  bins 

NZMA:  Number  of  horizontal  bins 

NCOSA:  Number  of  angular  bins 

NBATCH:  =1 

NHIS:  Number  of  neutrons  per  energy  band  run 

DELT:  Time  bin  widths 

VOLB:  Volume  of  the  spatial  bins 

The  remaining  fifteen  records  consist  of  D,  E,  F,  G. 

D  (dimensioned  360):  Sum  of  the  time  integrated  results 

E  (dimensioned  360):  Sum  of  the  squares  of  the  individual 
batch  time  integrated  results 
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F  (dimensioned  10400):  Sum  of  the  time  dependent  results 

G  (dimensioned  10400):  Sum  of  the  squares  of  the  individual 
batch  time  dependent  results 

The  results  are  stored  in  the  F  and  G  arrays  as  if  they  were 
dimensioned  (NTIMA,  NRADA,  NCOSA  +  NZMA).  Similarly  for  the 
D  and  E  arrays  (NRADA,  NCOSA  +  NZMA).  The  results  stored  in  the 
15  records  in  order  of  storage  is: 

Theta  and  vertical  currents  -  fast 

Radial  and  horizontal  currents  -  fast 

Ionization  (gammas)  -  fast 

Ionization  (neutrons)  -  fast 

Ionization  (neutrons  and  gammas)  -  fast 

Theta  and  vertical  currents  -  thermal 

Theta  and  vertical  currents  -  fast  and  thermal 

Radial  and  horizontal  currents  -  thermal 

Radial  and  horizontal  currents  -  fast  and  thermal 

Ionization  (gammas)  -  thermal 

Ionization  (gammas)  -  fast  and  thermal 

Ionization  (neutrons)  -  thermal 

Ionization  (neutrons)  -  fast  and  thermal 

Ionization  (neutrons  and  gammas)  -  thermal 

Ionization  (neutrons  and  gammas)  -  fast  and  thermal 


303 


I'iM  mm** 


{W? 1  •  rfi  * J«  '■  *J\  Ilf  J I  -1  »>  n™ 


REFERENCES 


1.  C.  A.  Stevens  and  R.  R.  Schaefer,  "Neutron  Induced 
Electrical  Sources,  "  AFWL-TR-69-8,  Air  Force  Weapons 
Laboratory,  Kirtland  AFB,  New  Mexico,  1968. 

2.  C.  A.  Stevens,  "Neutron  Induced  Electrical  Sources-II, " 
AFWL-TR-69-60,  Air  Force  Weapons  Laboratory, 

Kirtland  AFB,  New  Mexico,  1969. 

3.  C.  A.  Stevens  and  R.  E.  Dietz,  "Late-Time  EMP  Sources 
Including  Ground  Contributions,  "  AFWL-TR-71-8,  Air 
Force  Weapons  Laboratory,  Kirtland  AFM,  New  Mexico, 
1970. 

4.  E.  A.  Straker  and  M.  L.  Gritzner,  "Neutron  and  Secondary 
Gamma  Ray  Transport  in  Infinite  Homogeneous  Air, " 
ORNL-4464,  Oak  Ridge  National  Laboratory,  Oak  Ridge, 
Tennessee,  1969. 

5.  Private  Communication  from  V.  W.  Pine  of  AFWL. 

6.  Written  by  H.  M.  Murphy,  Jr.  of  AFWL. 

7.  H.  Kahn,  "Applications  of  Monte  Carlo,  "  RM-1237-AEC, 
RAND  Corporation,  Santa  Monica,  California,  1956. 

8.  H.  Goldstein,  Fundamental  Aspects  of  Reactor  Shielding, 
Addison -Wes ley  Publishing  Company,  Reading,  Mass.,  1959. 

9.  J.  H.  Ray,  G.  Grochowski,  and  E.  S.  Troubetzkoy, 

"Neutron  Cross  Sections  of  Nitrogen,  Oxygen,  Aluminum, 
Silicon,  Iron,  Deuterium,  and  Beryllium, "  TJNC-5139, 

United  Nuclear  Corporation,  White  Plains,  New  York,  1965. 

10.  M.  K.  Drake,  "Neutron  and  Gamma  Ray  Production  Cross 
Sections  for  Silicon,  "  DASA  2099,  GA-8628,  Gulf  General 
Atomic,  San  Diego,  California,  1968. 


304 


11.  S.  H.  Dupree,  H.  A.  Sandmeier,  G.  E.  Hansen,  W.  W. 
Engle,  Jr.,  and  F.  R.  Mynatt,  ’’Time-Dependent 
Neutron  and  Photon  Transport  Calculations  Using  the 
Method  of  Discrete  Ordinates,  "  LA  4557,  Los  Alamos 
Scientific  Laboratory,  Los  Alamos,  New  Mexico,  1970. 

12.  D.  C.  Irving,  R.  M.  Freestone,  Jr.,  and  F.  B.  K.  Kam, 
”05R,  A  General-Purpose  Monte  Carlo  Neutron  Transport 
Code,  "  ORNL-3622,  Oak  Ridge  National  Laboratory, 

Oak  Ridge,  Tennessee,  1965. 

13.  J.  Lindhard,  et  al. ,  ’’Integral  Equations  Governing 
Radiation  Effects  (Notes  on  Atomic  Collisions,  HI),  " 

Mat.  Fys.  Medd.  Dan.  Viol.  Selsk.  33,  No.  10,  1963. 


305 


